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Thrust Vector Control Heat Transfer Modeling 



Abstract 

The report presents heat transfer modeling of Thrust Vector control 
systems using the PHOENICS computer code* 

Simple two-dimensional wedge and blunt bodies have been examined in 
supersonic cold flow, for both- laminar and turbulent flow cases. 

The research presents a numerical solution of the supersonic compressible 
viscous two-dimensional flow field. Post calculations were done to estimate 
skin friction coefficient, surface heat flux, heat transfer coefficient and 
Stanton number distributions in both wedge and blunt cases. 
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Pressure 
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Heat flux 
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Gas constant [j/kg*k] 
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Time [S] 
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Temperature 
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Specific heat ratio 
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Boundary layer thickness 
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Dynamic viscosity [kg*m/s] 
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General exchange coefficient 
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Density [kg/m^ ] 




Constants used in turbulent model 




Any property at the grid node 
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1 , Introduction 



This report describes a numerical analysis of heat transfer of a typical 
jet vane configuration used for thrust vector control. The work was carried 
out under contract Nos. N62271-85-M-0443 and N6227 1--86-M-0206 , for the Naval 
Pos tgraduate School. 

The tasks to be accomplished under the first contract were: 

Task I: Formulate the conservation equations of momentum energy for 

two-dimensional, supersonic flow in geometries typical of thrust vector 
control systems. 

Task II : Formulate boundary conditions for these equations appropriate 

to thrust vector control systems. 

The tasks to be accomplished under the second contract were: 

Task I : Continue and update the formulation of thrust vector control 

geometries based on the input from the Naval Weapons Center 
(NWC). 

Task II : Construct the computational model for implementation in the 

PHOENICS code, of the thrust vector control geometries and 
flow conditions provided by NWC. 

Task III : Run the PHOENICS code for the previously formulated models. 

Analyze and interpret the PHOENICS results for surface 
temperature *and heat flux. 

Thrust vector control components such as jet vanes and jet tabs are 
exposed to high speed hot gases at the exit of a rocket nozzle. 

Estimation of the heat transfer from the hot exhaust gases to the vane is 
major consideration in the correct design of a vane, and its ability to 
survive during its mission. 
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The research work was done under the framework of the tasks. A brief 
survey of what has been done according to the task is given: 

Task 1 (M-0443): Heat transfer modeling of thrust vector control vane 

requires supersonic compressible viscous flow analysis* 

In order to meet the requirements, the conservation differential 
equations of mass momentum energy and the two k~e turbulent equations were 
formulated, and additional algebraic formulas for the relations between 
pressure density and the equation of state for ideal gas. 

Task II (M-0443): The physical dimensions of the flow field grid were 

chosen and the boundary conditions for the Navier-S tokes , energy and the two 
k-q turbulent model equations were given. 

Task I (M-0206): Working on the task, the actual configuration of a jet 

vane that is presently being tested at NWC has been modeled. The geometry 
being used is a wedge which has the same half angle and dimensions as the NWC 
jet vane. 

Task II (M-0206): BFC (Body fitted coordinate) version of PHOENICS code 

(Ref. 3) was used for calculating the flow-field and heat transfer over the 
model. Using the BFC, a better geometrical approximation to vane shape could 
be achieved. 

Non-Uniform grids have been utilized in order to model complicated regions 
in the flow field. Relaxation parameters and false times teps options were 
adjusted to enable efficient computer runs with good convergence. 

Task III (M-0206): In carrying out this task, four major runs have been 

analyzed : 

Two geometric configurations were used: wedge vane and blunt vane (see 

Figures 1, 2, 3, 4); each one in both laminar and turbulent flow conditions. ^ 

J 
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Numerical results for fluid field and thermodynamic properties of 
pressure, temperature, density, Mach number and velocities are given in 
appendix C. 

Post-calculations of heat transfer coefficient, skin friction coefficient 
and Stanton number are given in Figures (6, 7, 8, 9, 10, 11). 

The next chapters describe in more detail the process of building the 
model and the analysis of the results. 
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2. PHOENICS Description 



The present work addresses the heat transfer modeling of thrust vector 
control systems. In this effort the Navier-Stokes approach is applied by 
using a computer code which is capable of simulating a large number of fluid 
flow, heat transfer and chemical reaction processes which arise in industry 
and elsewhere. This code is cailed PHOENICS, which is an acronym standing 
for: ’Parabolic, Hyperbolic or Elliptic Numerical Integration Code Series.’ 

The name comes from the fact that the differential equations of fluid flow, 
etc. arise in forms classified by mathematicians as parabolic, hyperbolic or 
elliptic; and PHOENICS solves these equations, whatever their form. 

Built into PHOENICS are the major conservation laws of physics (mass, 
momentum, and energy) applied to a large number of continuous subdomains 
called ’cells,* into which the domain of study is artificially divided. The 
number of cells can be few or many according to the requirements of the 
problem. Because of numerical stability considerations the restrictions on 
cell refinement can become particularly burdensome in the calculation of a 
turbulent boundary layer where a very fine mesh near the wall may be 
required. 

When supplied with appropriate information concerning: the physical 

properties of the materials, the geometrical and other constraints, the inlet 
and/or initial conditions, PHOENICS computes the corresponding solutions to 
the relevant differential equations, expressing them as tables of numbers 
describing the field of velocity, temperature concentration, etc. 

Detailed information about PHOENICS is given in [Ref. 3]. 
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2*1 The Structural Principle of PHQENICS 



The code consists of three major parts: Satellite subroutine, 

Ground subroutine and Earth library. 

The satellite subroutine is the main input subroutine and should 
provide the answers to the questions: 

- what kind of process is to be simulated 

- what are the properties of the fluid 

- what are the shape and size of the domain 

- how fine is the grid to be employed 

- to what degree of accuracy is the calculation to be continued 

- and what output should be provided 

Ground subroutine is active during the computing process and is used for 
updating properties which vary with time, temperature, etc. For example: 
viscosity depends on temperature or density depends upon pressure and 
temperatures , etc. 

Earth library is the main solver generator. It is given as a binary 
library and does not enable the user access to the source code. 

2.2 Numerical Scheme 

The numerical scheme used by the code is the simpler (semi-implicit 
method for pressure-linked equations revised) (Ref. 9). The scheme was 
developed by Patankar, S. V. and Spalding, D. B. 

The scheme requires an additional dependent variable, the pressure 
correction, which has no physical meaning but should take part in the 
process . 

The value of the pressure correction should tend to zero in the 
convergence process . 

Two additional differential equations are solved: for the pressure, and 

for the pressure correction. 
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3 . Geometry and Dimensions 



Symmetrical 2-D planar geometry, which is shown in Figure 2, was chosen 
to be the approximation of the MWC vane in Figure 1. 

Two geometrical profiles were examined, one with wedge leading edge and 
the second with blunt leading edge. 

The dimensions of the domain in Figure 3 and 4 satisfy aspect ratio of 
10:1 in the vertical y coordinate. A high aspect ratio in the coordinate is 
important for the assumption of free stream conditions at the upper boundary. 
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4 . Assumptions 



Postulating the right or the wrong assumptions has the most influence on 
modeling process. The stage was carried out very carefully in order to make 
the most compatible model with reality. 

4.1 Steady state ; 

The modeling assumes steady state physical phenomenon process. 

- — (all properties) = 0 

O L 

This is a valid assumption since the time constant for the 
convection process is much shorter than the time constant for the wall 
conduction. 

By assuming the wall temperature to be constant, the two procedures 
are decoupled. 

In hot flow it is important to run the code for a wide range of wall 
temperature which will take into account the influence of different 
temperatures on the heat convection process. 

4.2 Cold Air Flow 

Ambient temperature air flow which was utilized by NWC experiments 
is being used in the computations. 

4.3 Ideal Gas 

The gas is assumed to satisfy the ideal gas equation of state 

p = pRT (4.1) 

This is a fairly good assumption for nonreactive gas flow. In spite of the 
values of static temperature can decrease to 200[k], the density remain 
relatively low. 
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This assumption is an important simplification to the solution in Ref. 10 
which used the isentropic relation between pressure and density instead 




4.4 Constant Pr,y : 

Prantdl number and y (ratio of specific heats) were found to have 
negligible variations in the temperature range of the model. (200k f 350k) 

4 .5 Varying Viscosity and Thermal Conductivity : 

y and k are much more dependent on temperature especially very close to 
the solid wall where values of y and k influence strongly the shear and heat 
transfer mechanism. To account for the temperature dependence power law 
relations have been formulated for y and k. 

'T 0.666 

M = MgCf — ) (A. 3) 

0.666 

k = k (-S — ) (4.4) 

^o 

4.6 Parallel Flow 

Gas flow at the exit of the exhaust nozzle is more likely to be a conic 
source flow than parallel flow. 

If the half angle of the nozzle is small, (a < 15°), parallel flow is a 
good assumption 

Negligible Radiation 

Assessments that were done showed that heat convection is at least one 
order of magnitude greater than heat flux by radiation. 
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4,8 Laminar and Turbulent Solutions 



In order to overcome lack of ability to predict transition, separated 
laminar and turbulent calculations were done for each case. The turbulent 
solution utilizes the (k-e) eddy viscosity model i^f, 5, 

4.9 Constant Wall Temperature 

The vane wall is assumed to have constant temperature during the time of 
calculation. 
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5 • Governing Equations 



The conservation equations for the compressible flow of the mathematical 
model consists of a viscous, Newtonian perfect gas consisting of the following 
six differential equations: 

Conservation of Mass: 

-^(p) + v(p^) = 0 (5.1) 

Conservation of momentum: 

-j|-Cp<f) + V Cp^(|) - uV4>) 7P (5.2) 

where 4) is V or W velocity component for y and z direction* 

Conservation of Energy 

(ph) + V(p^h - ^ Vh) = -gH (5.3) 

where h is the total enthalpy. 

h = Cp 

where To is the total temperature 
To = Tstat*(l 

In the case of laminar flow the governing equations (5.1), (5.2), (5.3) 
are sufficient to determine a solution when proper boundary conditions are 
applied and the equation of state (4.1) is provided. 

Turbulence Model: 

In turbulent flow it is necessary to hypothesize a turbulence idel 
relating the turbulent viscosity to the other problem variables. 



10 



The model used in PHOENICS is the eddy viscosity (k-g) model (Ref. 3, 

Ref. 5] . In this model k, the turbulent kinetic energy and e, the turbulence 
dissipation rate, are treated as properties of the flow and conservation 
equations are postulated for these properties. The two conservation equations 
are: one for k, the kinetic energy of turbulence: 



Dk _ 9 . ^eff 9k . 

Dt 9Xj '‘Ok 



+ Gk - e 



(5.4) 



Second equation for e, the dissipation rate of turbulence 

De ^ 9 . ^eff 9e 

Dt “ 9Xj'- og 3Xj 



D £ 9 / ff9£ \ , £/^|-i 

•5F “ sx7> * k =2^) 



(5.5) 



where 



Gw = vt (. 



9Ui 9Uj 9Ui 

) 



(5.6) 

(5.7) 



9Xj 9Xi' 9Xj 

‘'eff = ''lam ^ 

C 2 , a\/^f Og, are empirical constants which are provided in PHOENICS. 

The reason for using the (k-e) model is because it is the most verified 
model for engineering applications. It combines simplicity, universality, and 
realism of predictions in most cases. 

Two additional differential equations are solved also in order to satisfy 
the SIMPLER algorithm as was mentioned in chapter 2.2. The description of the 
pressure and pressure correction equations is provided by Ref. 9. 
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6 . Input Variables 



The properties of mach no, stagnation presence and temperature of the gas 
were provided by NWC; additional properties were taken from air tables: 

Mach number: 

= 3.2 

Stagnation pressure: 

Po = 55.1Q5 [Pa] 

Stagnation temperature : 

To = 555,55 [K] 

Gas constant 
R = 287. (J/kg*k] 

Specific heat ratio 
Y = 1.35 

Laminar Prandtl Number 
Pr => 0.7 

Turbulent Prandtl Number 
Pr^ = 0*9 

Constant Pressure Specific Heat 
Cp = R/(1-1 /y) (J/kg-k) 

Laminar Viscosity 
p => 0.1716 * 10~5 * (T/273.)0*^®® 

Thermal Conductivity 
k =• w Cp/Pf 

The gas properties in the inlet boundary are equivalent to the properties 
at nozzle exit. Inlet properties are calculated from the stagnation values in 
the combustion chamber. The calculation was done by assuming one dimensional 



12 



isentropic expansion from combustion chamber to the nozzle exit (inlet for 
the vane). 



Pressure 


Pi = Pq/CI 


(6.1) 


Temperature 


Tl = Tq/(1 + ^ m2) 


(6.2) 


Density 


Pi = P/Rt 


(6.3) 


Enthalpy 


hi = CpTi 


(6.4) 


Sonic Velocity 


Ci = /Trtt 


(6.5) 


Velocity 


Wi = C*M 


(6.6) 


subscript i signifies inlet property 
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7 . Boundary Conditions 



The flow field described in Figures 3, 4 has four boundaries, which can 
be named: inlet, outlet, frees tream boundary and solid wall. 

Super sonic flows have a hyperbolic mathematical nature. The field 
consists of influence zones, the flow at every point is governed only by its 
influence zone, basically by the upwind stream. 

As a consequence from the discussion, it’s obvious that the outlet 
boundary condition has no influence on the upstream flow. The boundary values 
that are given at the outlet are to satisfy some numerical needs only. 

7.1 Inlet 

Parallel uniform flow with known velocity, enthalpy, pressure, and 
density: equation (6.1), (6.3), (6.4), (6.6) are given at the left boundary 

of the grid. In PHOENICS this is specified as the LOW side of the first Z 
cell. 

In turbulent flow, boundary conditions are supplied for k and e. The 
values that are given are based on empirical values: 



where GH is half the vane thickness 
7 .2 Outlet 

As was mentioned previously, the outlet has negligible effect on the 
results. The only property that is specified at the outlet is the pressure. 



ki = 0.0 



(7.1) 



= 0.16 kl*5/(5*GH) 



(7.2) 
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7 .3 Frees tream Boundary 



Assuming that the upper boundary is chosen to be far enough away, the 
default boundary condition option of PHOENICS is used. This implies a line of 
symmpfry where all gradients are zero. 

7.4 Solid WalL 

Zero velocity and constant wall enthalpy (temperature) are assumed on the 
wall. In PHOENICS the wall is the SOUTH side of the first y cell. The high 
enthalpy and velocity gradients near the wall demands a refined grid close to 
the wall. Values of shear stress and heat flux are calculated to first order 
accuracy using: 



8W 



W, 



^ 3y “ ^ 



(7.1) 



qw 



ii_ 

Pr 3y Pr Ayi/ 2 



(7.2) 



In turbulent flow, a wall function is used to provide the wall condition 
for velocity, enthalpy, k, and e 



7 .5 Wall Function 

The wall problem in the numerical computation of flows, especially in 
turbulent flow, is an old one and most authors have adopted similar 
techniques. In effect they **6ridge over" the region very close to the wall by 
introducing special functions which are called wall functions. These are 
often empirical in origin. Accounts may be found in Ref. 11. 

The problem arises as follows. Turbulence dies out, close to the wall, 
because the no slip condition and the rigidity of the wall make all the 
velocity components fall to zero. The consequence is that the effective 
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viscosity and other transport properties fall there to their laminar values 
and the result is a rapid variation with distance from the wall both of the 
(j) ’s and of their gradients* 

Where signifies general dependent variable, It is possible to compute 
these variations in detail, by using a computer code such as PHOENICS on two 
conditions; 

(i) the grid points must be packed into the region of steep gradient 
changes closely enough for sufficient numerical accuracy to be obtained 

(ii) the functions appearing in the turbulence model equations must 
properly represent the influence of local Reynolds number on turbulence. 

Under the conditions above, the wall function sequences in the program 
act as follows: 

The Reynolds number is first evaluated, based on the resultant velocity 
parallel to the wall, on the distance from the wall to the grid node and on 
density and laminar viscosity. If this Reynolds number is less than 132.25 
(the value at which the laminar and turbulent wall function intersect) a 
laminar wall function is used. If this Reynolds number turns out to be 
greater than 132.25 the velocity variation is logarithmic and the 
corresponding shear stress coefficient is evaluated. This corresponds to the 
commonly used ’’log law’* wall function. [Ref. 4] 

7 .6 Boundary Conditions in Phoenics 

PHOENICS utilizes source terms for creating boundary conditions. The 
form of the source term of each dependent variable (j> is: 

= ([m] -h C^) (V^ - (^p) (7.3) 

where: m - is mass flux source 

(j) p - is the value of the dependent variable at point near the boundary 

C(J), V(j> - two coefficients specified by the user. The source term for 
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mass flux is simply 



Sm = (Vm - Pp) (7.4) 

where; Pp -is the pressure near the boundary and Cm* are two coefficients. 
The values of C(}i and V()i for the dependent variables in SATELLITE are: At the 

Inlet; 

(7.5) 

Vm = Po Pi/^o (7.6) 

Cw = Ch = Ck = Cg = 0. (7.7) 

Vw = Wi (7.8) 

Vh = hi (7.9) 

Vk = Ki (7.10) 

Ve = £i (7.11) 

At the Outlet: 

Cni = 1000 * Wi*pi/Pi (7.12) 

Vm = Pi (7.13) 

At the Wall (laminar) 

Cw = u/(0.5Aui) (7.14) 

Vw = o (7.15) 

Ch = u/Pr/(0.5*Aui) (7.16) 

Vh = V^w (7.17) 

At the Wall (turbulent) 

= Ck ™ — WALL (7.18) 

V„ = Vk = Vg = 0 (7.19) 

Vh = Cp*Tw 
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(7.20) 



8. Mesh Generation 



In this work a two-dimensional mesh is being used with 18 x 29 cells in 
the y and z coordinate respectively* A Nonuniform grid has been used for 
both directions. Figures 3 and 4 shows the grid in the z direction. A finer 

grid is used in the blunt region, IZ = (7 x 17), and in the zone, where the 

inclined wall transitions to a straight wall, IZ =» (23 r 26). 

In the y coordinate, except in the boundary layer region, the grid is 
uniform. To obtain a finer grid resolution in the boundary layer for the 
laminar flow case the first five cells in the y direction from the wall obey 
the following proportionality relationship: 

TV 3 Amfl'v 

BYFRAC (lY) = (r-j—) (8.1) 

Where BYFEIAC(IY) is the distance from the south side to the north side of the 

cell of particular interest, divided by total length of the domain, lY is the 

cell number, A^ax maximum allowable cell height, and GH is the half 
thickness of the TVC jet vane. 

A fine grid resolution for the turbulent flow case is set up in the same 
way as laminar flow. The only difference comes from the selection of the 
first five cells in y direction. The following calculation shows the 
difference . 

From the laminar solution and the given properties the following are 
known: 

w =» 885.2[m/s] 

Wlam 1*10”5 [N.s/m] 

Po = 5.5*106[Pa] 

Pstatic = 1.048*105(Pa) 

Y = 1.35 
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p = 1.835 [kg/m] 



Using the values above and the length of vane, which is 0.095m, A 

corresponding Reynolds number was calculated: 

_ P« w« Z (1.835*888.5*0.095) _ , * .,^^6 

• 

Using a power law correlation for the boundary layer thickness: 

- = 0.37 * Rez"l/5 



( 8 . 2 ) 



From equation (8.2) the boundary layer thickness at the high end of the 
domain has been calculated as 5 » 2 * 10“^ [m] 

Ay 

With Re based on W„ the velocity parallel to the wall, — ^ the distance 

2 

from the wall to the first grid node, p„ the density, and the laminar 

viscosity. Ay must satisfy the condition 

poo '^oo 



^A 



> 132.25 



or 



2 Ulam 

Therefore the interval of Ay is chosen such that 
2 * 10"3(m] > AY > 6.48 * 10"^(m] 



AY > 6.48 * 10"6[m] 



In this effort using the relationship 
BYFRAC(IY) . (-^)' (^) 

Ay has been calculated as Ay = 8 * 10”^ [m] which is in the required interval. 

For both the laminar and turbulent cases, cells in the z direction were 
adjusted so that, the points where possible physical phenomena such as shock 
waves and expansion fans are expected, very fine cells were used. In the 
other parts of the domain larger cells were used. 
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9 • Heat Transfer Analysis 



Skin friction and heat transfer quantities were calculated in both 
laminar and turbulent cases and they are shown in Figures (6 t 11). 

9 • 1 Laminar Calculation 

In laminar flow fluxes can. be derived directly from the gradients near 
the wall* The first cell is close "enough” to the wall and gradients of 
velocity and enthalpy do not change much in this region near the wall* The 
shear stress and heat flux in the laminar case will be: 






(7.1) 



_ \i ^1 - 

AY772 



(7.2) 



The skin friction coefficient and Stanton number will be: 




(9.1) 



St = <lw/ 1 P»lJco( hf - h„)] 
where h^. is the recovery enthalpy 



(9.2) 



hr 1 + 
ho ■ 1 + ( y-1 ) 




(9.3) 






2 



00 



r- is the recovery factor 



r ■ /Pr (laminar flow) 



(9.4) 
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The coefficient of heat transfer in convection was calculated using 



P« Cp St 



(9.5) 



9.2 Turbulent Calculations 

In turbulent flow the gradients of velocity and enthalpy near the wall 
are very steep and change rapidly with distance from the wall. 

Direct calculation of flux gradients is not accurate in this case. The 
log law approach is used to calculate skin friction. In the calculations 
using PHOENICS flow field, the following relation has been used. 



C - 2 Pw 
^ w2p®3.33 

oo 

To obtain equation 
a starting point. 



(9.6) 

9.6, the turbulent kinetic energy equation has been used as 
[Ref. 5], 



Dk 3 -Vk 

Dt 8y 8y' 



(9.7) 



, , r bt ^3 u 2 k 1 

+ k ( 1— (^TT-) - — J 



k '3y ' ut 

The source term of the turbulent kinetic energy equation should be zero near 
the wall which means 



_ P^ k 

k ^ 3y ^w 






= 0 



(9.8) 



therefore the shear stress on the wall can be defined as: 



Tw = Pw ’^w ^9.9) 

where is the turbulent kinetic energy on the wall, is the density on the 
wall and Cj) = 0*09 [Ref* 5], substituting the values above into the Blasius 
skin friction relation the Cf equation becomes: 
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2 



(9.10) 



Cf = 



poo 




00 




oo 



The heat transfer quantities are evaluated from the Chilton-Colburn form of 
Reynolds analogy. 

St = (Cf/2) * Pr‘2/3 (9.11) 

q„ = St*p„*U„*(hr-h^) (9.12) 

where equation (9.3) is used to evaluate hj. with the recovery factor given as: 
1/3 

r = Pj- (turbulent flow) (9.13) 

The convective heat transfer coefficient is calculated by using equation (9.5) 
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10. Code and Computer 



PHOENICS 81, Body Fitted Coordinate (BFC) version has been used in the 
computations (see Ref. 3). PHOENICS has been installed on NPS IBM 3033 MVS 
1.3 computer. 400 sweeps per computer run provided a reasonable convergence 
in all runs except the turbulent blunt case continuity error of less than 
4.10”** has been achieved in thd three runs. 

The continuity error is the total summation of the absolute mass 
imbalance in all cells divided by the inlet mass flux. CPU time consumption 
varies from case to case as follows: 



Laminar Wedge 


630 


CPU 


Seconds 


Turbulent Wedge 


630 


CPU 


Seconds 


Laminar Blunt 


630 


CPU 


Seconds 


Turbulent Blunt 


1542 


CPU 


Seconds for 1000 sweeps 
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1 1 . Results and Discussion 



The results of the calculations are available on appendix c. The tabular 
results include the values of pressure, velocities, enthalpy, temperature mach 
number, density, turbulent kinetic energy and rate of turbulent dissipation. 
The values are given in 18 x 29 cells points. 

Skin friction and heat transfer results are shown in Figures (5-11). 
Laminar and turbulent skin friction and Stanton number in wedge flow show 
improvement compared to the results reported by Yukselen (Ref. 10). The lines 
are smoother and the oscillations at the end were eliminated. Basically the 
magnitudes are similar to those in Ref. 10. 

Laminar blunt values are similar except near the beginning. The 
beginning, as expected in blunt zone, creates higher rates of heat transfer. 
Even though the blunt geometry used is a multi-wedge shape it should predict 
the correct values except for the stagnation point itself. 

Turbulent blunt skin friction has different behavior. It has a very 
large value at the first point and then undershoots to values that are smaller 
than for wedge. It should also be kept in mind that the convergence of this 
case wasn't very successful. 
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12. Conclusions and Recoaimendations 



1. PHOENICS was found to be a friendly code for simulating complicated 
mixed heat transfer fluid dynamics problems, 

2. Derivation of heat transfer properties to a vane solid wall in 
laminar and turbulent flow has been installed in the code. It can be used for 
predictions of heat transfer rate in both cold and hot gas flow. 

3. Two features have been added to the code in NFS: The restart option 

and the use of initial field, make it possible to simulate time dependent 
processes and solve the temperature variation in the vane itself. 
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Figure 1= NWC Jet vane configuration. 
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Figure 2 - NWC Jet Vane Approximation 
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Figure 3= Wedge vane domain and grid 
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Figure 4- Blunt vone domain and grid. 
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Figure 5: No. along the Vane 
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Figure 6= in Laminar flow. 
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Figure 7 - in Turbulent flow 
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Figure 9= Turbulent stanton number. 
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Figure 11* Coefficient of heat convection in turbulent flow. 
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Appendix A 



Satellite Listing 

Two subroutines SATELLITE and GROUND had to be changed and improved. The 
full list is enclosed in Appendix A and B. 

VAN4SAT and VANTSAT are the laminar and turbulent SATELLITE subroutines, 
the first has the blunt geometry and the second has the wedge (it can be 
changed easily from wedge to blunt and vice versa) VAN4GRD and VANTGRD are 
exactly the same. They are the GROUND subroutines, VANTGRD is given in 



Appendix B 



FILEi vantsat fortran A1 



C$DIRECTIVExxSATLIT AMI LEITNER 

C LAMINAR SOLUTION FOR NHC5 NY=18 NZ=29 YN=GTH 

C LECSAT CONVERTED TO DIAMSAT 

C XFILE NAME: MODBFCST.FTN 

C ^ABSTRACT: SATELLITE MODEL MAIN PROGRAM. THIS VERSION IS 
C FOR USE WITH THE BODY-FITTED COORDINATE SCHEME (SUMMER 1984 
C VERSION) PROVIDED AS AN ATTACHMENT TO SPRING 1983 PHOENICS. 

C XDOCUMENTATION: PHOENICS INSTRUCTION MANUAL (SPRING 1983) 

C WITH BODY-FITTED COORDINATES INSTRUCTION SUPPLEMENT 

C (SUMMER 1984). 

C ^AUXILIARY SUBROUTINES (TAPES, ETC.) ARE IN SATELLITE LIBRARY 

C SERVICEU, WHICH MUST BE INCLUDED IN LINK EDIT TO RUN. 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 1 STARTS: 





CHAPTER 1 COMMON BLOCKS AND USER'S DATA. 

C 

INCLUDE (CMNGUS) 

INCLUDE (CMNGRF) 

INCLUDE (GUSSEQ) 

COMMON/CPI/IPWRIT, IDUM(243) 

DIMENSION GDTAPE(3),DFAULT(4) 

DIMENSION ARRAY1(309),ARRAY2(194),ARRAY3(421) 

LOGICAL ARRAYl , LSPDA, WRT, RD, NAMLST 
INTEGER ARRAY2,XPLANE,YPLANE,ZPLANE 

INTEGER P1,PP,U1,U2,V1,V2,W1,W2,R1,R2,RS,EP,H1,H2,H3,C1,C2, 
&C3,C4 

REAL NORTH, LOW 
LOGICAL BFC 

EQUIVALENCE (ARRAYl ( 1 ), CARTES) , (ARRAY2(1 ) , NX) 

EQUIVALENCE ( ARRAY3( 1 ) , SPAREl ( 1 ) ) , ( Ml , R1 ) , (M2 , R2) 

EQUIVALENCE ( LSTRUN, INTGR( 12) ),( NAMLST, L0GIC(88) ) 

EQUIVALENCE ( LOGIC! 20 ), BFC) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 1 ENDS. 
C$DIRECTIVEXXCMNBF1$$ 

C THIS FILE CONTAINS SATELLITE COMMON BLOCKS FOR BFC'S 
C FI MUST BE DIMENSIONED TO GREATER THAN OR EQUAL TO 
C (NX+NY+17XNZ+24xNXxNY+6x(NX+1)x(NY+1)+6XND) . THE VALUE 
C OF THE DIMENSION MUST BE SET AS NBFC IN GROUP 6 OF SATLIT. 
COMMON/F0B/FK5000) 

COMMON/CIB/ND/CIC/KOORD 

COMMON/CID/KDBGG, KDBGMF, KDBGCD, KDBIND, KDBMFX, KDBCDT, KDBPCS, 

8 KDBGUV,KDBGPV 

COMMON/CI E/KDBGS , KDBI NS 
COMMON/CIF/IGEN/CIG/NCART 

C THE FOLLOWING ARRAYS MUST BE EXACTLY DIMENSIONED FOR NXPl, 

C NYPl AND NZPl, BUT MAY BE OVER DIMENSIONED FOR ND. 

C THE BFRAC ARRAYS MUST BE DIMENSIONED TO ALLOW FOR SETTINGS 
C IN SATLIT, THEY MAY BE OVER DIMENSIONED. 

COMMON/CRA/XW( 19 , 30 , 1 )/CRB/XE( 19,30,1) 

& /CRC/YS(2,30,1)/CRD/YN(2,30,1) 

& /CRE/ZL(2,19,1)/CRF/ZH(2,19,1) 

& /CRG/RCON/CRH/DARCY/CRI/BXFRAC( 99 )/CRJ/BYFRAC( 99 ) 

& /CRK/BZFRAC(99) 

C0MM0N/CLA/ST0RSA(6 ) , STORWD(6 ) , STORP, STORPE, STORPN, 

& ST0RPH,ST0R1,ST0R2,ST0R3,ST0UNV,PRTBFC,ST0CRN 

COMMON/CLC/BFPLOT 

LOGICAL STORP, STORPE, STORPN, STORPH, STORl , ST0R2, ST0R3, 

& STORSA, STORWD, STOUNV, PRTBFC, BFPLOT, STOCRN 

C END 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 1 STARTS: 

C GRAFFIC ARRAYS DIMENSIONED AS NEEDED... 

COMMON/GRAFl/PHIl ( 1 ) /GRAF2/PHI2( 1 ) 

C . POROSITY 8 SPECIAL DATA ARRAYS DIMENSIONED AS NEEDED... 
DIMENSION PE(1,1,1),PN(1,1,1),PH(1,1,1),PC( 1,1,1) 

DIMENSION LSPDA( 1 ) , ISPDA( 1 ) , RSPDA( 1 ) 

C USER PLACES HIS VARIABLES, ARRAYS, EQUIVALENCES ETC. HERE. 
EQUIVALENCE(RAIR,RE(21)),(GAMA,RE(22)),(GSWP,RE(23)) 
1,(GPR,RE(24)),(TW,RE(25)),(GEMU1,RE(26)),(JEMU1,INTGR(D) 

C USER PLACES HIS DATA STATEMENTS HERE. 

DATA NLSP,NISP,NRSP/1,1,1/ 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 1 ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 2 STARTS: 



VANOOOlO 

VAN00020 

VAN00030 

VAN00040 

VAN00050 

VAN00060 

VAN00070 

VAN00080 

VAN00090 

VANOOlOO 

VANOOllO 

VAN00120 

VAN00130 

VAN00140 

VAN00150 

VAN00160 

VAN00170 

VAN00180 

VAN00190 

VAN00200 

VAN00210 

VAN00220 

VAN00230 

VAN00240 

VAN00250 

VAN00260 

VAN00270 

VAN00280 

VAN00290 

VAN00300 

VAN00310 

VAN00320 

VAN00330 

VAN00340 

VAN00350 

VAN00360 

VAN00370 

VAN00380 

VAN00390 

VAN00400 

VAN00410 

VAN00420 

VAN00430 

VAN00440 

VAN00450 

VAN00460 

VAN00470 

VAN00480 

VAN00490 

VAN00500 

VAN00510 

VAN00520 

VAN00530 

VAN00540 

VAN00550 

VAN00560 

VAN00570 

VAN00580 

VAN00590 

VAN00600 

VAN00610 

VAN00620 

VAN00630 

VAN00640 

VAN00650 

VAN00660 

VAN00670 

VAN00680 

VAN00690 

VAN00700 

VAN00710 

VAN00720 
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FILE* VANTSAT FORTRAN A1 



C 

CHAPTER 2 SET CONSTANTS. AND ARRANGE FILE MANIPULATIONS. 

C 

C PLEASE DO NOT ALTER. OR RE-SET. ANY OF THE REMAINING 
C STATEMENTS OF THIS CHAPTER. 

DATA CELL. EAST.WEST.NORTH. SOUTH. HIGH. LOW. VOLUME/ 

& 0..!.. 2. .3. .A. .5. .6. .7. / 

DATA P1.PP.U1.U2.V1.V2.W1.W2.R1.R2.RS.KE.EP.H1.H2.H3.C1.C2. 
SC3 . CA/ 1^2. 3. A. 5. 6. 7. 8. 9. 10. 11. 12. 13. lA. 15. 16. 17. 18. 19. 2 0/ 

DATA FIXFLU. FIXVAL.ONLYMS. WALL/1. E-1 0.1 .El 0.0.0. -10.0/ 

DATA IPLANE.XPLANE.YPLANE.ZPLANE/0. 1.2.3/ 

DATA WRT. RD. DFAULT/ . TRUE. . . FALSE . . AHDEFA. 4HULT .;. ^HDTA/ . IHG/ 
DATA GDTAPE/AHGUSI.AHE1.D.2HTA/ 

DATA NLDATA.NIDATA.NRDATA/309.194.A21/ 

DATA NLCREG.NTCVRG/60.350/ 

DATA TITPP.TITCl.TITC2.TITC3/3HRH0.<iHMACH.AHTEMP. AHCFST/ 

CALL TAPESC lO.GDTAPE. 3. 1 . A)<NRDATA) 

c read default file if blockdata absent 

IF(INTGR1(29).NE.10) GO TO 2 

CALL WRIT40CA0HDATA ESTABLISHED IN BLOCK DATA. ) 

GO TO 3 

2 CALL DEFLT 

CD 2 CALL TAPESC1.DFAULT.A.2.4XNRDATA) 

CD CALL DATAIO(RD.l) 

CALL WRITAOCAOHDATA TAKEN FROM DEFAULT. DTA ON GROUP A/C) 

3 CALL WRITAOCAOHFILE MODSTL.FTN IS THE SATLIT USED. ) 

L0GIC(89) = .TR'JE. 

C 

CHAPTER 3 DEFINE DATA FOR NRUN RUNS. 

C 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 2 ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 2 STARTS* 

C— GROUP AIMULTI-RUNS * RUN( 1-30X .T . . 29X . F. > 

C 

RUN(1)=. FALSE. 

C NOTE* ALL RUNS ARE DEACTIVATED AT THIS POINT - USER SHOULD 
C ==== SWITCH ON ONE ONLY OF RUNS 1-A IN NEXT STATEMENT. 
RUN(1)=.TRUE. 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 2 ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 3 STARTS* 
DO 10 IRUN=1.30 

IF( .NOT.RUN(IRUN)) GO TO 10 
NRUN=NRUN+1 
LSTRUN=IRUN 
10 CONTINUE 

DO 999 IRUN=1.LSTRUN 

IF( .NOT.RUN(IRUN)) GO TO 999 
INTGR(ll) = IRUN 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 3 ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 3 STARTS: 

C— ALL INTEGER VARIABLES ARE DEFAULTED TO 0. AND REAL VARIABLES 
C TO 0.0. UNLESS OTHERWISE INDICATED. 

C E.G. BY VARIABLE<10>. OR <10. 0> AS APPROPRIATE. 

C THE DEFAULT SETTINGS OF ALL LOGICAL VARIABLES ARE ALWAYS 

C INDICATED. E.G. VARIABLE< .T . >. OR VARIABLE< . F . > . 

C 

C— RUNl 

C 

C— GROUP 1. FLOW TYPE * 

C PARAB<.F.>.CARTES<.T.>.ONEPHS<.T.> 

C 

C— GROUP 2. TRANSIENCE * 

C STEADY<.T.>.ATIME.LSTEP<1>.FSTEP<1> 

C ' TLAST<l.E10>.TFRAC(l-30X30)a.> 

C SERVICE SUBROUTINE FOR 'NT* POWER-LAW TIME STEPS* 

C CALL GRDPWRCO. NT. TLAST. POWER) 

C 

C— - GROUP 3. X-DIRECTION * 

C NX<1>.XULAST<1 .0>.XFRAC(l-30) 

C SERVICE SUBROUTINE FOR POWER-LAW GRID* 

C CALL GRDPWRCl. NX. XULAST. POWER) 

C 



VAN00730 

VAN007A0 

VAN00750 

VAN00760 

VAN00770 

VAN00780 

VAN00790 

VAN00800 

VAN00810 

VAN00820 

VAN00830 

VAN008A0 

VAN00850 

VAN00860 

VAN00870 

VAN00880 

VAN00890 

VAN00900 

VAN00910 

VAN00920 

VAN00930 

VAN009A0 

VAN00950 

VAN00960 

VAN00970 

VAN00980 

VAN00990 

VANOIOOO 

VANOlOlO 

VAN01020 

VAN01030 

VANOIOAO 

VAN01050 

VAN01060 

VAN01070 

VAN01080 

VAN01090 

VANOllOO 

VANOlllO 

VAN01120 

VAN01130 

VANOllAO 

VAN01150 

VAN01160 

VAN01170 

VAN01180 

VAN01190 

VAN01200 

VAN01210 

VAN01220 

VAN01230 

VAN01240 

VAN01250 

VAN01260 

VAN01270 

VAN01280 

VAN01290 

VAN01300 

VAN01310 

VAN01320 

VAN01330 

VAN013A0 

VAN01350 

VAN01360 

VAN01370 

VAN01380 

VAN01390 

VANOIAOO 

VANOIAIO 

VAN01A20 

VAN01A30 

VANOIA^O 
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ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo ooooo oooo 



FILE: VANTSAT FORTRAN A1 



group 4. Y-DIRECTION : 

NY<1>/YVLAST<1 . 0> , YFRAC( 1-30) , RINNER, SNALFA 
SERVICE SUBROUTINE FOR POWER-LAW GRID: 

CALL GRDPWR(2, NY, YVLAST, POWER) 

NY=18 



GROUP 5. Z-DIRECTION : 

NZ<1>,ZWLAST<1 . 0>,ZFRAC(1-30) 

SERVICE SUBROUTINE FOR POWER-LAW GRID: 
CALL GRDPWR(3,NZ,ZWLAST, POWER) 

NZ=29 



— GROUP 6. MOVING GRID OR DISTORTED (BODY-FITTED-) GRID : 
— MOVING GRID : 

MGRID,IZW1,IZW2,AZW2,BZW2,CZW2,PINT,ZW2M1T 



— BODY-FITTED GRlD — 
BFC<.T.>,IGEN<1>,ND<1>,NBFC<5000>,KOORD,RCON 
BXFRACd-NXXl .0,NXM1*0.0> 

BYFRAC(1-NYX1.0,NYM1X0.0> 

BZFRAC(1-NZX1.0,NZM1*0.0> 

SERVICE SUBROUTINE FOR SUB-DOMAIN SPECIFICATION (FOR IGEN=1 
ONLY) : 

CALL DOMAIN(ID,IXF,IXL,IYF, IYL,IZF, IZL) 
XE(l-NYPl,l-NZPl,l-NDX(NYPlXNZPlXND)3a.0>, 

XW( 1-NYPl , 1-NZPl , 1-ND) , 

YN(1-NXP1,1-NZP1,1-NDX(NXP1XNZP1*ND)*1.0>, 

YS(1-NXP1,1-NZP1,1-ND), 

ZH( 1-NXPl , 1-NYPl, 1-NDX(NXP1KNYP1XND)X1 . 0>, 

ZL ( 1-NXPl , 1-NYPl , 1-ND) ,5TORSA( 1-6 X6X . F .>, STORWD( 1-6 )<6* . F . >, 
STORP<.F.>,STORPE<.F.>,STORPN<.F.>,STORPH<.F.>,STOUNV<.F.>, 
PRTBFC< . F .>, DARCY, BFPLOT< . F . > 

CYCLIC BOUNDARY CONDITIONS ARE DEFAULTED INACTIVE ; 

TO ACTIVATE THEM AT SELECTED IZ SLABS USE SERVICE SUBROUTINE: 
CALL XCYIZdZ, .TRUE. ) 

SERVICE SUBROUTINE TO DEACTIVATE CURVATURE TERMS IN U, V 
AND W EQUATIONS ASSOCIATED WITH CURVATURE OF IX, lY, IZ 
GRID LINES RESPECTIVELY: 

CALL UCURVEdZ, .FALSE.) 

CALL VCURVEdZ, .FALSE.) 

CALL WCURVEdZ, .FALSE. ) 

NCART<1> 

XWARNINGSI I I I I I 

A) WHEN USING BFCS ST0VAR(H3), ST0VAR(C4), ST0VAR(21) ARE 
AVAILABLE ONLY FOR STORING NON-ORTHOGONAL VELOCITY 
COMPONENTS. 

B) MULTI-RUNS ARE NOT ALLOWED WITH BFC OPTION. 

C) MOVING GRID, TWO-PHASE AND PARABOLIC OPTIONS ARE NOT 
AVAILABLE WITH BFC OPTION. 

D) KE-EP TURBULENCE MODEL SHOULD BE USED WITH BFCS ONLY 
WHEN THE MAIN FLOW IS IN THE IZ DIRECTION. 

E) BUILT-IN GRAVITY TERMS DO NOT TAKE ACCOUNT OF BFCS. 
XNOTES 



A) THE STANDARD VELOCITY-FIELD PRINTOUT FOR THE 
VELOCITY RESOLUTES IS ACTIVATED IN THE USUAL 
WAY. AN ADDITIONAL OPTION EXISTS FOR PRINTING THE 
CARTESIAN VELOCITY-COMPONENTS WHICH MAY BE 
ACTIVATED BY SETTING THE FOLLOWING LOGICALS: 

ST0VAR(U2)=.T. FOR U-COMPONENT (CARTESIAN) 
ST0VAR(V2)=.T. FOR V-COMPONENT (CARTESIAN) 
ST0VAR(W2)=.T. FOR W-COMPONENT (CARTESIAN) 

SIMILARLY PRINTOUT OF NON-ORTHOGONAL VELOCITY 
COMPONENTS MAY BE ACTIVATED AS FOLLOWS: 

ST0VAR(C4)=.T. FOR U-COMPONENT (NON-ORTHOG) 
ST0VAR(H3)=.T. FOR V-COMPONENT (NON-ORTHOG) 
ST0VAR(21)=.T. FOR W-COMPONENT (NON-ORTHOG) 

B) BFC (TO ACTIVATE THE BFC OPTION), IGEN (THE CODE FOR METHOD 

OF GRID SPECIFICATION), ND (NUMBER OF SUB-DOMAINS) AND 
NBFC (THE FI ARRAY DIMENSION), MUST BE SET BEFORE 
"STANDARD BFC SECTION 2". ====== 
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C ALL OTHER BFC DATA MUST BE SET AFTER "STANDARD BFC 

C SECTION 2. ===== 

C C) NXPl, NYPl, NZPl STORE NX+1, NY+1, NZ+1; THESE ARE 

C AVAILABLE TO USER AFTER STANDARD BFC SECTION 2. 

C D) FOR IGEN=1 USE BXFRAC, BYFRAC & BZFRAC IN PLACE OF 

C XFRAC^YFRAC & ZFRAC. 

C 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD BFC SECTION 1 STARTS: 
C DEFAULT SETTINGS: 

NCART=10 

BFC=,TRUE. 

IGEN=1 

ND=1 

NBFC=5000 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD BFC SECTION 1 ENDS.: 

C XUSER SETS BFC, IGEN, ND AND.NBFC HERE: 

C 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD BFC SECTION 2 STARTS: 
CALL SBAI ( NXPl , NX+1 , NYPl , NY+1 , NZPl , NZ+1 ,1,0) 

IF(BFC) CALL BFCDFTC NBFC, XE, XW, YN, YS, ZH, ZL , ND, NXPl , NYPl , 
a NZP1,NZ) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD BFC SECTION 2 ENDS. 

C XUSER SETS ALL OTHER BFC VARIABLES HERE: 

C FUSING NONIFORM GRID 1-8 
GTH=65.E-3 
GTL=150.E-3 
GBETA=A. 

GBETA=GBETA)^3.1A15927/180 

GTAB=TAN(GBETA) 

DELMAX=2.E-3 

GNBL=5. 

GPWR=2. 

DO 6A IY=1,5 

64 BYFRAC(IY) = (FLOAT(IY)/GNBL)X)^GPWR)^DELMAX/GTH 
BYFRAC(6)=BYFRAC(5)+3.E-3/GTH 

DEL=(1 .-BYFRAC(6))/(FL0AT(NY)-GNBL-1) 

DO 65 IY=7,NY 

65 BYFRAC(IY)=BYFRAC(IY-1)+DEL 

BZFRAC(l)=10.E-3 
DO 66 IZ=2,5 

66 BZFRAC(IZ)=10 . E-3+BZFRACC IZ-1 ) 

BZFRAC(6)=BZFRAC(5)+5.E-3 

DO 67 IZ=7,9 

67 BZFRAC(IZ)=BZFRAC(IZ-l)+2.E-3 
DO 68 IZ=10,10 

68 BZFRAC(IZ)=BZFRAC(IZ-1)+1 .E-3 
DO 77 IZ=11,14 

77 BZFRAC(IZ)=BZFRAC(IZ-l)+.5E-3 
DO 78 IZ=15,15 

78 BZFRAC(IZ)=BZFRAC(IZ-1)+1 .E-3 
BZFRACC 16 ) =BZFRAC( 1 5 ) +2 . E-3 
BZFRACC 17 ) =BZFRAC( 16 )+3 . E-3 
BZFRACC 1 8 ) =BZFRAC( 1 7 )+5 . E-3 
DO 69 IZ=19,22 

69 BZFRACC IZ)=BZFRAC( IZ-1 )+10. E-3 
BZFRACC 23 )=BZFRACC 22 )+3. E-3 
BZFRACC 24 )=BZFRACC 23 )+2. E-3 
BZFRACC 25 )=BZFRACC 24 )+2. E-3 
BZFRACC26) =BZFRACC25)+3 . E-3 
BZFRACC 27 ) =BZFRACC 26 ) +5 . E-3 

DO 71 IZ=28,NZ 

71 BZFRACC IZ)=BZFRACC IZ-1 )+10 . E-3 
DO 72 IZ=1,NZ 

72 BZFRACC IZ)=BZFRACCIZ)/GTL 
CALL D0MAINC1,1,NX,1,NY,1,NZ) 

DO 61 IX=1,NXP1 

DO 62 IY=1,NYP1 
ZLCIX,IY,1)=0.0 
62 ZHCIX,IY,1)=GTL 
DO 63 IZ=1,NZP1 
YNCIX,IZ,1)=GTH 
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63 YS(IX,IZ,1)=0.0 

c YSdX/13,1) SHOULD COME AFTER 
DO 662 IZ=16,25 

^662 YS(IXdZa) = (BZFRAC(IZ-l)-BZFRAC(3))xGTABXGTL 
DO 663 IZ=13,15 

GZ12=(BZFRAC(IZ-1)-BZFRAC(11))XGTL 
663 YS(IX,IZ,1)=SQRT(YS(IX,16,1)xGZ12x2.-GZ12xx2) 
DO 664 IZ=26,NZ 
664 YS(IX,IZ,l)=YSaX,25,l) 

61 CONTINUE 

ST0RSA(IFIX( LOW) )=. TRUE. 
STORSA(IFIX(HIGH))=.TRUE. 
ST0RSA(IFIX(S0UTH))=.TRUE. 
STORWD(IFIX(SOUTH))=.TRUE. 

STORP=.TRUE. 

PRTBFC=.TRUE. 

CDAR DARCY=1.E10 



GROUP 7. BLOCKAGE: BL0CK< . F. >, IPLANE, IPWRIT 

XSET CONSTANT POROSITIES OVER SUB-DOMAINS USING: 

CALL CONPORdR, TYPE, VALUE, IXF,IXL,IYF,IYL,IZF,IZL), WHERE: 
IR=RUN SECTION NUMBER, E.G. 1 FOR RUNl SECTION; *TYPE’= EAST, 
WEST, NORTH, SOUTH, HIGH, LOW & CELL. 'VALUE’ =WANTED POROSITY 
OVER REGION IXF, . . .IZL. 

XDIMENSION ARRAYS PE(NX, NY, NZ) , PN( NX, NY, NZ) , PH( NX, NY, NZ) , & 
PC(NX,NY,NZ) ABOVE. 

XFOR FULLY-BLOCKED CELLS (IE. 'VALUE'* 0.0) USER NEED SET ONLY 
THE 'CELL' POROSITY (TO ZERO), AS CELL-FACE AREAS ARE THEN 
AUTOMATICALLY ZEROED. 

XFOR SATELLITE PRINTOUT OF ALL POROSITIES IN DOMAIN, 'IPLANE'* 
XPLANE YPLANE OR ZPLANE, FOR DESIRED CROSS-SECTION DIRECTION. 

XFOR EACH 'TYPE' A MAXIMUM OF 10 CALLS TO CONPOR IS ALLOWED, 
BUT IF REQUIREMENTS EXCEED THIS PROVISION SET BLOCK*. T. & 
IPWRIT*-!, AND SET POROSITY ARRAYS EXPLICITLY HERE AS WANTED. 
IN THIS CASE, THE USER MUST SET A L L ELEMENTS OF 
ARRAYS PE, PN, PH, PC (MANY MAY BE 0.0 OR 1.0). HE MAY USE: 
CALL CR(PARRAY,VALUE,IXF,IXL,IYF,IYL,IZF,IZL,NX,NY,NZ) 

ANY NUMBER OF TIMES, TO SET 'PARRAY' (= PE, ETC.) TO 
'VALUE' OVER RANGE IXF TO IXL, lYF TO lYL, IZF TO IZL. 

XCONPOR MUST NOT BE USED IN CONJUNCTION WITH EXPLICIT 
SETTINGS OF THE ARRAYS (INCLUDING SETTINGS VIA CR). 



— GROUP 8. DEPENDENT VARIABLES TO BE SOLVED FOR OR STORED : 

S0LVAR(1-25)<25X.F.>,ST0VAR( 1-25 )<25X.F.>,C0NC1( 1-4 )<4X.T.> 
USE FOLLOWING NAMED INTEGERS FOR ARRAY ELEMENTS 1-20: 
P1,PP,U1,U2,V1,V2,W1,W2,M1,M2,RS,KE,EP,H1,H2,H3,C1,C2,C3,C4. 
S0LVAR(P1)*.TRUE. 

SOLVAR(PP)=.TRUE. 

S0LVAR(V1)*.TRUE. 

S0LVAR(W1)*.TRUE. 

S0LVAR(H1)*.TRUE. 

SOLVAR(KE)*.TRUE, 

SOLVAR(EP)*.TRUE. 

ST0VAR(V2)*.TRUE. 

ST0VAR(W2)*.TRUE, 

ST0VAR(C1)*.TRUE. 

ST0VAR(C2)*.TRUE. 

ST0VAR(C3)*.TRUE. 



— GROUP 9. VARIABLE LABELS : 

TITLEd-25X2HPl,2HPP,2HUl,2HU2,2HVl,2HV2,2HWl,2HW2,2HRl, 
2HR2, 2HRS, 2HKE, 2HEP, 2HH1 , 2HH2, 2HH3, 2HC1 , 2HC2, 
2HC3,2HC4,2HRX,2HRY,2HRZ, 2x4Hxxxx> 

TITLE(C1)*TITC1 

TITLE(C2)*TITC2 

TITLE(C3)*TITC3 

TITLE(PP)*TITPP 



GROUP 10 PROPERTIES: 

IRH01<1>,IRH02<1>,RHOK1.0>,RH02<1.0>, 
ARHOKl . 0>, BRHOKl . 0>,CRHOK1 . 0> 
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IEMUK1>. EMUKl . 0>, EMULAM<1 . E-10> 
IHSAT,H1SAT,H2SAT,PSATEX<1.0> 

SIGMA(l-25Xl.0,2.0a. ,l.E10,l..l .ElOa.U.ElO, 
Axl. 0,1.31^, 1.0,1 . ElO, 10X1. 0> 

IRH01=-1 

PT0T=55.E5 

T0T=555.55 

RAIR=287. 

GAMA=1.35 

CP=RAIR/(1-1/GAMA) 

TW=323. 

HWALL=TWXCP 

htot=cpxtot 

RHT0T=PT0T/T0T/RAIR 

L0GIC(87)=.TRUE. 

ARH01=RHT0T/PT0TXX( 1/GAMA) ^ 

BRH01=1./GAMA 
TURBULENT OR LAMINAR 
IEMU1=2 
IEMU1=-1 
JEMU1=IEMU1 
EMUl=l.E-5 
EMULAM=EMU1 
GEMU1=EMU1 
GPR=.7 

SIGMA(2A)=GPR 

SIGMA(14)=.9 



— GROUP 11 INTER-PHASE TRANSFER PROCESSES i 

ICFIP,CFIPS,IMD0T,CMD0T,CA1K1 .E6>,CA2K1.E6> 



GROUP 12 SPECIAL SOURCES : 

ISPCSOC 1-25), AGRAVX,AGRAVY,AGRAVZ,ABUOY, HREF 



— - GROUP 13 INITIAL FIELDS : 
FIINITC 1-25X25x1 .E-10> 

MACH NO. OF FREE STREAM 
GMACH=3.2 

A=1+(GAMA-1)/2XGMACHXX2 

TE=TOT/A 

RHE=RHT0T/AXX(1/(GAMA-D) 

PSTAT = PTOT/Axx(GAMA/(GAMA-D) 

RH01=ARH01XPSTATXXBRH01 

SONIC=SQRT(GAMAxRAIRxTE) 

WIN=S0NICXGMACH 

RKEIN=0.01xWINXx2 

EPIN=0 . 16XRKEINXX1 . 5/GTH/2 . 

FIINIT(W1)=WIN 

FIINIT(P1)=PSTAT 

FIINIT(H1)=HT0T 

FIINIT(KE)=RKEIN 

FIINIT(EP)=EPIN 



— GROUP lA BOUNDARY/INTERNAL CONDITIONS : 

IL00P1,IL00PN,XCYCLE<.F.>,PBAR,REGI0N(1-10X10X.T.> 

XN.B. ALL 10 REGIONS ARE DEFAULTED .TRUE.. THE USER SHOULD 
SET REGION(I)=. FALSE. FOR UNUSED REGIONS •!*. 

DO lA 1=1,10 
1^ REGION(I)=. FALSE. 



— - GROUP 15 TO 24; REGIONS 1 TO 10 

— ONLY THOSE REGIONS ARE ACTIVE WHICH ARE SPECIFIED BY THE 
. USER, PREFERABLY BY WAY OFi- 
CALL PLACE(IREGN,TYPE,IXF,IXL,IYF,IYL,IZF,IZL) 8 
CALL COVAL(IREGN,VARBLE,COEFF, VALUE) 

CALL PLACE(l,LOW,l,NX,l,NY,l,l) 

C CALL C0VAL(1,M1,FIXFLU,WINXRHE) 

CDAR CALL COVAL ( 1 , Ml , 1 . E-20 , 1 . E+20XWINXRHE) 
GCM=2XGAMA/WIN/(GAMA-1) 

GVM=PT0TXRHE/RHT0T 
CALL C0VAL(1,M1,GCM,GVM) 

CALL COVAL(l,Wl,ONLYMS,WIN) 
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c 



c 



CALL COVAL(1»H1,ONLYMS,HTOT) 

CALL C0VAL(1,KE,0NLYMS,RKEIN) 

CALL C0VAL(1,EP»0NLYMS,EPIN) 

CALL PLACE(2,HIGH,1,NX,1»NY,NZ,NZ) 

CALL COVAL(2,M1,FIXVAL,PSTATXO. ) 

CALL COVAL(2,M1»1000xWINxRHE/PSTAT,PSTAT) 
CALL C0VAL(2,H1,0NLYMS,HT0T) 

WALL ALONG THE VANE IZ(11,NZ) 

GCM=EMUl/( .5XBYFRAC(1)*GTH) 
DY1=BYFRAC(1)XGTH 
GOEFF=EMU1/(0.5 xDY1) 
GOEFH=EMU1/(0.5XDY1*SIGMA(2A)) 

CALL PLACE(3,S0UTH,1,NX,1,1,12,NZ) 

CALL COVAL(3,W1,GOEFF,0. ) 

CALL C0VAL(3,H1,G0EFH,HWALL) 

CALL C0VAL(3,W1,WALL, 0 . ) 

CALL C0VAL(3,H1,WALL,HWALLX- 
CALL COVAL(3,KE,WALL,0. ) 

CALL C0VAL(3,EP,WALL,0. ) 



— GROUP 25 GROUND STATION « 

GR0STA< . F . >, NAMLST< . F , > 

XNAMLST ACTIVATES NAMELIST IN GROUND. 
GROSTA=.TRUE. 



— GROUP 26 SOLUTION TYPE AND RELATED PARAMETERS « 
WHOL EP< . F . > , SUBPST< , F . >, D0NACC< . F . > 
WHOLEP=.TRUE. 



— GROUP 27 SWEEP AND ITERATION NUMBERS « 

FSWEEP<1>,LSWEEP<1>,LITHYD<1>,LITC<1>,LITKE<1>,LITH<1>, 
LITER( 1-25X9x1, -1,15X1> 

IVELF<1>,NVEL<1>,IVELL<10000>, 

IKEF<1>,NKE<1>,IKEL<10000>, 

IENTF<1>,NENT<1>,IENTL<10000>, 

ICNCF<1>,NCNC<1>,ICNCL<10000>, 

IRH01F<1>,NRHOK1>,IRH01L<10 00 0>, 

IRH02F<1>,NRH02<1>,IRH02L<10000> 

LSWEEP=1201 

GSWP=LSWEEP 

FSWEEP=801 

LITER(PP)=20 

LITER(V1)=5 

LITER(W1)=5 

LITHYD=2 



— GROUP 28 TERMINATION CRITERIA : 

ENDIT( 1-25X9X1. E-1 0,0. 5, 15X1. E-10> 
ENDIT(l)=l.E-5 



— GROUP 29 RELAXATION : 

RLXP<1 . >, RLXPXY<1 . >, RLXPZ<1 . >, RLXRH0<1 . >, RLXMDT<1 . >, 
DTFALS( 3-25X23X1. E10> 

DTFALS(Wl)=l.E-5 

DTFALS(Vl)=l.E-5 

DTFALS(KE)=l.E-5 

DTFALS(EP)=l.E-6 

RLXP=.3 



— GROUP 30 LIMITS i 

VELMAX<1.E10>,VELMIN<-1.E10>,RHOMAX<1.E10>,RHOMIN<1.E-10>, 
TKEMAX<1 . E10>,TKEMIN<1 . E-10>, EMUMAX<1 . E10>, EMUMIN<1 . E-10>, 
EPSMAX<1.E10>,EPSMIN<1.E-10>,AMDTMX<1.E10>,AMDTMN<-1.E10> 
EPSMAX=1.E13 



— GROUP 31 SLOWING DEVICES « SL0RH0<1 . >, SL0EMU<1 . > 
SL0RH0=.2 



— GROUP 32 PRINT-OUT OF VARIABLES « 

PRINT(1-25X.T., .F.,23X.T.>,SUBWGR<.F.> 
PRINT(C1)*.TRUE. 

PRINT(C2)=.TRUE. 
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PRINT(C3)=.TRUE. VAN05050 

PRINT(PP)=.TRUE. VAN05060 

C VAN05070 

C GROUP 33 MONITOR PRINT-OUT : VAN05080 

C IXMON<l>,IYMON<l>» IZMON<l>,NPRMON<l>,NPRMNT<l> VAN05090 

NPRMON=10 VAN05100 

IYM0N=2 VAN05110 

I2M0N=12 VAN05120 

Q VAN05130 

C GROUP 3A FIELD PRINT-OUT CONTROL i VAN051A0 

C NPRINT<100>,NTPRIN<100>,NXPRIN<1>,NYPRIN<1>,NZPRIN<1>, VAN05150 

C I2PRF<1>,ISTPRF<1>,IZPRL<10000>,ISTPRL<10000> - VAN05160 

C NUMCLS<10>,KOUTPT VAN05170 

NPRINT=LSWEEP VAN05180 

C VAN05190 

C— GROUP 35 TABLE CONTROL i VAN0520D ■ 

C TABLES<.F.>,NTABLE,NTABVR,LINTAB,NPRTAB,NMON, VAN05210 

C ITAB(l-8),MTABVR(l-8) VAN05220 

C VAN05230 

C GROUP 36-38 ARE NOT DOCUMENTED IN THE INSTRUCTION VAN052<i0 

C MANUAL AND ARE INTENDED FOR MAINTENANCE PURPOSES ONLY VAN05250 

C— GROUP 36 DEBUG PRINT-OUT SLAB AND TIME-STEP i VAN05260 

C IZPR1<1>,IZPR2<1>,ISTPRK1>,ISTPR2<1> VAN05270 

C VAN05280 

C— GROUP 37 DEBUG SWEEP AND SUBROUTINES : VAN05290 

C KEMU,KMAIN,KINDEX,KGEOM,KINPUT,KSODAT,KCOMPF,KSORCE, VAN05300 

C KS0LV1,KS0LV2,KS0LV3,KC0MPP,KADJST,KFLUX,KSHIFT,KDIF, VAN05310 

C KCOMPU,KCOMPV,KCOMPW,KCOMPR,KWALL,KDBRHO<-l>,KDBEXP,KDBMDT VAN05320 

C KDBGEN VAN05330 

C VAN05340 

C— GROUP 38 MONITOR, TEST, AND FLAG : VAN05350 

C M0NITR<.F.>,FLAG<.F.>,TEST<.T.>,KFLAG<1> VAN05360 

C END OF MAINTENANCE-ONLY SECTION VAN05370 

C VAN05380 

C GROUP 39 ERROR AND RESIDUAL PRINT-OUT « VAN05390 

C IERRP<1000>,RESREF(1,3-2A)<25X1.>,RESMAP<.F.>, VAN05A00 

C RESID(1-25X2X.F. ,23X.T.>,K0UTPT VAN05A10 

RESREF(1)=WINXRHE VAN05A20 

RESREF(7)=WINXRESREF(1) VAN05A30 

RESREF(5)=WINxrESREF(1)X0.1 VAN05AA0 

RESREF(H1)=HT0TXRESREF(1) VAN05A50 

RESREF(KE)=RKEINXRESREF(1) VAN05A60 

RESREF(EP)=EPINXRESREF(1) VAN05470 

IERRP=LSWEEP/20 VAN05480 

KOUTPT=LSWEEP/20 VAN05490 

C VAN05500 

C— - GROUP 40 SPECIAL DATA : LOGICC 1 . . 10 ) , INTGR( 1 . . 10 ) , RE( 21 . . 30 ) , VAN05510 

C NLSP<1>,NISP<1>,NRSP<1>,SPDATA<.F.>,LSPDA(1),ISPDA(1),RSPDA(1) VAN05520 

C USE FIRST 10 ELEMENTS OF ARRAYS LOGIC & INTGR AND 21ST VAN05530 

C TO 30TH OF ARRAY RE FOR TRANSFERRING SPECIAL DATA FROM VAN05540 

C SATELLITE TO GROUND, BUT IF REQUIREMENTS EXCEED THIS VAN05550 

C PROVISION SET SPDATA = .T., AND DIMENSION ARRAYS LSPDA, VAN05560 

C ISPDA, RSPDA ABOVE AND IN GROUND AS NEEDED, AND SET HERE VAN05570 

C VAN05580 

C— GROUP 42 RESTARTS AND DUMPS : SAVEM< . F. >, RESTRT< . F. >, KINPUT VAN05590 

SAVEM=.TRUE. VAN05600 

BFPLOT=.TRUE. • VAN05610 

RESTRT=.TRUE. VAN05620 

C VAN05630 

C VAN05640 

C— GROUP 43 GRAFFIC i VAN05650 

C GRAPHS<.F.>,0RTH0G<.T.>,ANTSYM,NPRT<1>,ITITL<5X4HXXXX> VAN05660 

C— FOR A GRAFFIC RUN, DIMENSION PHIl S PHI2 AS FOLLOWSi VAN05670 

C PHIKNXXNYXNZXNM) VAN05680 

C PHI2((NX+2)X(NY+2)X(NZ+2)X(NM+IBLK)) , WHERE VAN05690 

C NM=NO. OF VARIABLES STORED + DENSITY(-IES) VAN05700 

C IBLK=0 IF BLOCK=.FALSE.,=4 IF A 3D RUN, VAN05710 

C =3 IF A 2D.YZ RUN. VAN05720 

C VAN05730- 

IF(IRUN.EQ.l) GO TO 900 VAN05740 

900 CONTINUE VAN05750 

C ALL RUNS VAN05760 
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C- 

C 



C 

C 

C 



CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 3 ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION A STARTS. 

WRITEGENERAL DATA ON TO THE GUSIEl.DTA TAPE, ETC... 

IF(SPDATA) CALL WRTSPC( LSPDA, NLSP, ISPDA, NISP, RSPDA, NRSP) 
IF(BLOCK) CALL WRTPORCPE, PN, PH, PC, NX, NY, NZ, IPLANE) 

IF(BFC) CALL WRTBFC(1A,NBFC,XE,XW,YN,YS,ZH,ZL, 
aND,NX+l,NY+l,NZ+l,NZ,PRTBFC) 

OLD PRACTICES RETAINED FOR REFERENCE. 

IFCSPDATA) CALL SPCDAT(IRUN) 

IF(BLOCK) CALL PORDAT(IRUN) 

IF(GRAPHS) CALL SORT(IRUN) 

IF(RESTRT) GO TO 902 
DO 901 INDVAR=1,25 

IF(IFIX(FIINIT(INDVAR)+0.1) .NE. 10101) GO TO 901 
CALL FLDDAT(IRUN) 

GO TO 902 

901 CONTINUE 

902 CALL DATAI0(WRT,10) 

IF(MONITR) CALL DATAIO( WRT, -6 ) 

999 CONTINUE 
STOP 
END 

CXXX IGEN=1 SO BFCXYZ NOT REQUIRED. 

C*K* COMMENT OUT BOTH VERSIONS. 



SUBROUTINE BFCXYZ (NXPl , NYPl , NZPl ) 

RETURN 

END 
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C$DIRECTIVEXXSATLIT AMI LEITNER 

C LAMINAR SOLUTION FOR NWC5 NY=18 NZ=29 YN=GTH 

C LECSAT CONVERTED TO DIAMSAT 

C xFILE NAME: MODBFCST.FTN 

C ^ABSTRACT: SATELLITE MODEL MAIN PROGRAM. THIS VERSION IS 
C FOR USE WITH THE BODY-FITTED COORDINATE SCHEME (SUMMER 1984 
C VERSION) PROVIDED AS AN ATTACHMENT TO SPRING 1983 PHOENICS. 

C ^DOCUMENTATION: PHOENICS INSTRUCTION MANUAL (SPRING 1983) 

C WITH BODY-FITTED COORDINATES INSTRUCTION SUPPLEMENT 

C (SUMMER 1984) . 

C ^AUXILIARY SUBROUTINES (TAPES, ETC.) ARE IN SATELLITE LIBRARY 

C SERVICEU, WHICH MUST BE INCLUDED IN LINK EDIT TO RUN. 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 1 STARTS: 



C 

CHAPTER 1 COMMON BLOCKS AND USER'S DATA. 

C 

INCLUDE (CMNGUS) 

INCLUDE (CMNGRF) 

INCLUDE (GUSSEQ) 

C0MM0N/CPI/IPWRIT,IDUM(243) 

DIMENSION GDTAPE(3),DFAULT(4) 

DIMENSION ARRAYK309), ARRAY2(194),ARRAY3(421) 

LOGICAL ARRAYl , LSPDA, WRT, RD, NAMLST 
INTEGER ARRAY2,XPLANE,YPLANE,ZPLANE 

INTEGER P1,PP,U1,U2,V1,V2,W1,W2,R1,R2,RS, EP, HI , H2, H3, Cl , C2, 
8C3,C4 

REAL NORTH, LOW 
LOGICAL BFC 

EQUIVALENCE ( ARRAYl ( 1 ), CARTES) , (ARRAY2( 1 ), NX) 

EQUIVALENCE ( ARRAY3( 1 ) , SPAREl ( 1 ) ) , (Ml , R1 ) , (M2, R2) 

EQUIVALENCE ( LSTRUN, INTGR( 12) ) , ( NAMLST, L0GIC(88)) 

EQUIVALENCE ( LOGIC( 20 ) , BFC) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 1 ENDS. 
C$DIRECTIVEXXCMNBF1$$ 

C THIS FILE CONTAINS SATELLITE COMMON BLOCKS FOR BFC'S 
C FI MUST BE DIMENSIONED TO GREATER THAN OR EQUAL TO 
C (NX+NY+17xNZ+24xnXxNY+6x(NX+1)x(NY+1)+6xND) . THE VALUE 
C OF THE DIMENSION MUST BE SET AS NBFC IN GROUP 6 OF SATLIT. 
COMMON/F0B/FK5000) 

COMMON/CIB/ND/CIC/KOORD 

COMMON/CID/KDBGG,KDBGMF,KDBGCD,KDBIND,KDBMFX,KDBCDT,KDBPCS, 

Sl KDBGUV,KDBGPV 

COMMON/CIE/KDBGS,KDBINS 
COMMON/CIF/IGEN/CIG/NCART 

C THE FOLLOWING ARRAYS MUST BE EXACTLY DIMENSIONED FOR NXPl, 

C NYPl AND NZPl, BUT MAY BE OVER DIMENSIONED FOR ND. 

C THE BFRAC ARRAYS MUST BE DIMENSIONED TO ALLOW FOR SETTINGS 
C IN SATLIT, THEY MAY BE OVER DIMENSIONED. 

C0MM0N/CRA/XW(19, 30, 1 )/CRB/XE( 19 , 30 , 1 ) 

& /CRC/YS(2,30, 1)/CRD/YN(2,30,1) 

& /CRE/ZL(2,19,1)/CRF/ZH(2,19,1) 

& /CRG/RC0N/CRH/DARCY/CRI/BXFRAC(99)/CRJ/BYFRAC(99) 

8 /CRK/BZFRAC(99) 

COMMON/CLA/STORSA ( 6 ) , STORWD( 6 ) , STORP , STORPE, STORPN , 

& STORPH, STORl , ST0R2, ST0R3, STOUNV, PRTBFC, STOCRN 

COMMON/CLC/BFPLOT 

LOGICAL STORP, STORPE, STORPN, STORPH, STORl , ST0R2, ST0R3, 

& STORSA,STORWD, STOUNV, PRTBFC, BFPLOT, STOCRN 

C END 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 1 STARTS: 

C GRAFFIC ARRAYS DIMENSIONED AS NEEDED... 

COMMON/GRAFl/PHIl ( 1 ) /GRAF2/PHI2( 1 ) 

C POROSITY & SPECIAL DATA ARRAYS DIMENSIONED AS NEEDED... 
DIMENSION PE(1,1,1),PN(1,1,1),PH(1,1,1),PC( 1,1,1) 

DIMENSION LSPDA(1),ISPDA(1),RSPDA(1) 

C USER PLACES HIS VARIABLES, ARRAYS, EQUIVALENCES ETC. HERE. 
EQUIVALENCE(RAIR,RE(2D), (GAMA,RE(22)), (GSWP,RE(23)) 

1, (GPR,RE(24)), (TW,RE(25)), (GEMU1,RE(26)), (JEMU1,INTGR(D) 

C USER PLACES HIS DATA STATEMENTS HERE. 

DATA NLSP,NISP,NRSP/1,1,1/ 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 1 ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 2 STARTS: 
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CHAPTER 2 SET CONSTANTS, AND ARRANGE FILE MANIPULATIONS, 



C- 

C 

C 



CD 

CD 



please do not ALTER, OR RE-SET, ANY OF THE REMAINING 
STATEMENTS OF THIS CHAPTER. 

DATA CELL, EAST, WEST, NORTH, SOUTH, HIGH, LOW, VOLUME/ 
a 0.,1.,2.,3.,A.,5.,6.,7. / 

DATA P1,PP,U1,U2,V1,V2,W1,W2,R1,R2,RS,KE,EP,H1,H2,H3,C1,C2, 
aC3,C4/l, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 15, 14,15,16, 17, 18, 19, 20/ 

DATA FIXFLU, FIXVAL, ONLYMS, WALL/1. E-10,1. £10,0.0,-10.0/ 

DATA IPLANE,XPLANE,YPLANE,ZPLANE/0,1,2,3/ 

DATA WRT,RD,DFAULT/.TRUE., . FALSE ., 4HDEFA, 4HULT . >4HDTA/, IHG/ 
DATA GDTAPE/4HGUSI,4HE1.D,2HTA/ 

DATA NLDATA, NI DATA, NRDATA/309, 194, 421/ 

DATA NLCREG,NTCVRG/60,350/ 

DATA TITPP,TITC1,TITC2,TITC3/3HRH0,4HMACH,4HTEMP,4HCFST/ 

CALL TAPES(10,GDTAPE,3, 1,4XNRDATA) 

read default file if blockdata absent 

IF(INTGR1(29) .NE.IO) GO TO 2 

CALL WRIT40(40HDATA ESTABLISHED IN BLOCK DATA. ) 

GO TO 3 
CALL DEFLT 

CALL TAPES(1,DFAULT,4,2,4XNRDATA) 

CALL DATAI0(RD,1) 

CALL WRIT40(40HDATA TAKEN FROM DEFAULT. DTA ON GROUP A/C) 

CALL WRIT40(40HFILE MODSTL.FTN IS THE SATLIT USED. ) 

L0GIC(89)=.TRUE. 



C 

CHAPTER 3 DEFINE DATA FOR NRUN RUNS. 

C 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 2 ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 2 STARTS: 

C— GROUP 41MULTI-RUNS s RUN( 1-30X . T . , 29X. F . > 

C 

RUN(1)=. FALSE. 

C NOTE: ALL RUNS ARE DEACTIVATED AT THIS POINT - USER SHOULD 
C ==== SWITCH ON ONE ONLY OF RUNS 1-4 IN NEXT STATEMENT. 
RUN(1)=.TRUE. 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 2 ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 3 STARTS: 
DO 10 IRUN=1,30 

IF( .NOT.RUN(IRUN)) GO TO 10 
NRUN=NRUN+1 
LSTRUN=IRUN 
10 CONTINUE 

DO 999 IRUN=1,LSTRUN 

IF( .NOT.RUN(IRUN)) GO TO 999 
INTGR(ll) = IRUN 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 3 ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 3 STARTS: 

C— ALL INTEGER VARIABLES ARE DEFAULTED TO 0, AND REAL VARIABLES 
C TO 0.0, UNLESS OTHERWISE INDICATED. 

C E.G. BY VARIABLE<10>, OR <10. 0> AS APPROPRIATE. 

C THE DEFAULT SETTINGS OF ALL LOGICAL VARIABLES ARE ALWAYS 
C INDICATED, E.G. VARIABL E< . T . >, OR VARIABL E< . F . > . 

C 

C RUNl 

C 



C GROUP 1. FLOW TYPE : 

C PARAB<.F.>,CARTES<.T.>,ONEPHS<.T.> 

C 

C— GROUP 2. TRANSIENCE : 

C . STEADY<.T.>,ATIME,LSTEP<1>,FSTEP<1> 

C TLAST<1.E10>,TFRAC(1-30X30X1.> 

C SERVICE SUBROUTINE FOR 'NT* POWER-LAW TIME STEPS: 

C CALL GRDPWR(0, NT, TLAST, POWER) 

C 

C GROUP 3. X-DIRECTION : 

C NX<l>,XULAST<1.0>,XFRAC(l-30) 

C SERVICE SUBROUTINE FOR POWER-LAW GRID: 

C CALL GRDPWRd, NX, XULAST, POWER) 

C 
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VAN4SAT FORTRAN A1 



GROUP Y-DIRECTION : 

NY<1>,YVLAST<1 .0>,YFRAC(1-30),RINNER,SNALFA 
SERVICE SUBROUTINE FOR POWER-LAW GRIDi 
CALL GRDPWRC2. NY, YVLAST, POWER) 

NY=18 



— GROUP 5. Z-DIRECTION : 

NZ<1>,ZWLAST<1 .0>,ZFRAC(l-30) 

SERVICE SUBROUTINE FOR POWER-LAW GRIDi 
CALL GRDPWR(3,NZ,ZWLAST, POWER) 

NZ=29 



— GROUP 6. MOVING GRID OR DISTORTED (BODY-FITTED) GRID * 
MOVING GRID : 

MGRID,IZW1,IZW2,AZW2,BZW2,CZW2,PINT,ZW2M1T 



— BODY-FITTED GRID — 
BFC<.T.>,IGEN<1>,ND<1>,NBFC<5000>,KOORD,RCON 
BXFRACC l-NXXl . 0, NXMIXO . 0> 

BYFRACd-NYXl .0,NYM1X0 .0> 

BZFRACd-NZXl .0,NZM1X0.0> 

SERVICE SUBROUTINE FOR SUB-DOMAIN SPECIFICATION (FOR IGEN=1 
ONLY) s 

CALL DOMAINdDdXF, IXL , lYF, lYL , IZF, IZL ) 
XE(l-NYPld-NZPl,l-NDX(NYPl3(NZPlXND)3a .0>, 

XWd-NYPld-NZPl, 1-ND), 

YN(1-NXP1,1-NZP1,1-NDX(NXP1)^NZP1XND)X1 .0>, 
YS(1-NXP1,1-NZP1,1-ND), 

ZH(l-NXPld-NYPl,l-NDX(NXPlXNYPlXND)^^l .0>, 

ZL(1-NXP1,1-NYP1,1-ND),ST0RSA(1-6)<6X.F.>,ST0RWD(1-6X6X.F.>, 
STORP<.F.>,STORPE<.F.>,STORPN<.F.>,STORPH<.F.>,STOUNV<.F.>, 
PRTBFC<.F.>, DARCY, BFPL0T<.F.> 

CYCLIC BOUNDARY CONDITIONS ARE DEFAULTED INACTIVE ; 

TO ACTIVATE THEM AT SELECTED IZ SLABS USE SERVICE SUBROUTINE* 
CALL XCYIZdZ, .TRUE. ) 

SERVICE SUBROUTINE TO DEACTIVATE CURVATURE TERMS IN U, V 
AND W EQUATIONS ASSOCIATED WITH CURVATURE OF IX, lY, IZ 
GRID LINES RESPECTIVELY: 

CALL UCURVEdZ, .FALSE. ) 

CALL VCURVEdZ, .FALSE.) 

CALL WCURVEdZ, .FALSE. ) 

NCART<1> 

^WARNINGSI I II I I 

A) WHEN USING BFCS ST0VAR(H3), STOVAR(C^), ST0VAR(21) ARE 
AVAILABLE ONLY FOR STORING NON-ORTHOGONAL VELOCITY 
COMPONENTS. 

B) MULTI-RUNS ARE NOT ALLOWED WITH BFC OPTION. 

C) MOVING GRID, TWO-PHASE AND PARABOLIC OPTIONS ARE NOT 
AVAILABLE WITH BFC OPTION. 

D) KE-EP TURBULENCE MODEL SHOULD BE USED WITH BFC’S ONLY 
WHEN THE MAIN FLOW IS IN THE IZ DIRECTION. 

E) BUILT-IN GRAVITY TERMS DO NOT TAKE ACCOUNT OF BFCS. 
)^NOTES 



A) THE STANDARD VELOCITY^FI EL D PRINTOUT FOR THE 
VELOCITY RESOLUTES IS ACTIVATED JN THE USUAL 
WAY. AN ADDITIONAL OPTION EXISTS FOR PRINTING THE 
CARTESIAN VELOCITY-COMPONENTS WHICH MAY BE 
ACTIVATED BY SETTING THE FOLLOWING LOGICALSi 

ST0VAR(U2)=.T. FOR U-COMPONENT (CARTESIAN) 
ST0VAR(V2)=.T. FOR V-COMPONENT (CARTESIAN) 
ST0VAR(W2)=.T. FOR W-COMPONENT (CARTESIAN) 

SIMILARLY PRINTOUT OF NON-ORTHOGONAL VELOCITY 
COMPONENTS MAY BE ACTIVATED AS FOLLOWS: 

STOVAR(C^)=.T. FOR U-COMPONENT (NON-ORTHOG) 
ST0VAR(H3)=.T. FOR V-COMPONENT (NON-ORTHOG) 
ST0VAR(21)=.T. FOR W-COMPONENT (NON-ORTHOG) 

B) BFC (TO ACTIVATE THE BFC OPTION), IGEN (THE CODE FOR METHOD 

OF GRID SPECIFICATION), ND (NUMBER OF SUB-DOMAINS) AND 
NBFC (THE FI ARRAY DIMENSION), MUST BE SET BEFORE 
"STANDARD BFC SECTION 2". ====== 
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C ALL OTHER BFC DATA MUST BE SET AFTER "STANDARD BFC 

C SECTION 2. ===== 

C C) NXPl, NYPl, NZPI STORE NX+1, NY+1, NZ+1; THESE ARE 

C AVAILABLE TO USER AFTER STANDARD BFC SECTION 2. 

C D) FOR IGEN=1 USE BXFRAC, BYFRAC & BZFRAC IN PLACE OF 

C XFRAC,YFRAC & ZFRAC. 

CXXXXXXXXXXXXXX^^ STANDARD BFC SECTION 1 STARTS: 

C DEFAULT SETTINGS: 

NCART=10 

BFC=.TRUE. 

IGEN=1 

ND=1 

NBFC=5000 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD BFC SECTION 1 ENDS.: 

C XUSER SETS BFC, IGEN, ND AND NBFC HERE: 

C 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD BFC SECTION 2 STARTS: 
CALL SBAI (NXPl, NX+1, NYPl, NY+1, NZPI, NZ+1, 1,0) 

IF(BFC) CALL BFCDFT(NBFC, XE, XW, YN, YS, ZH, ZL, ND, NXPl, NYPl , 
a NZP1,NZ) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD BFC SECTION 2 ENDS. 

C XUSER SETS ALL OTHER BFC VARIABLES HERE: 

C XUSING NONIFORM GRID 1-8 
GTH=65.E-3 
GTL=150.E-3 
GBETA=4. 

GBETA=GBETAX3. 1415927/180 
GTAB=TAN(GBETA) 

DELMAX=2.E-3 

GNBL=5. 

GPWR=4. 

DO 64 IY=1,5 

64 BYFRAC(IY)=(FLOAT(IY)/GNBL)XXGPWRXDELMAX/GTH 
BYFRAC(6)=BYFRAC(5)+3.E-3/GTH 
DEL=(1.-BYFRAC(6))/(FL0AT(NY)-GNBL-1) 

DO 65 IY=7,NY 

65 BYFRAC(IY)=BYFRAC(IY-1)+DEL 

C ZZ 

BZFRAC(l)=10.E-3 
DO 66 IZ=2,5 

66 BZFRAC(IZ)=10.E-3+BZFRAC(IZ-l) 

BZFRAC( 6 ) =BZFRAC( 5)+5 . E-3 

DO 67 IZ=7,9 

67 BZFRAC(IZ)=BZFRAC(IZ-l)+2.E-3 
DO 68 IZ=10,11 

68 BZFRAC(IZ)=BZPRAC(IZ-l)+.5E-3 
BZFRAC(12)=BZFRAC(ll)+l.E-3 
DO 77 IZ=13,14 

77 BZFRAC(IZ)=BZFRAC(IZ-l)+.5E-3 
DO 78 IZ=15,15 

78 BZFRAC(IZ)=BZFRAC(IZ-1)+1 .E-3 
BZFRAC(16)=BZFRAC(15)+l.E-3 
BZFRAC( 17 )=BZFRAC( 16 )+2 . E-3 
BZFRAC(18)=BZFRAC(17)+7 .E-3 
DO 69 IZ=19,22 

69 BZFRAC(IZ)=BZFRAC(IZ-l)+10.E-3 
BZFRAC(23)=BZFRAC(22)+3.E-3 
BZFRAC(24)=BZFRAC(23)+2.E-3 
BZFRAC(25)=BZFRAC(24)+2.E-3 
BZFRAC(26)=BZFRAC(25)+3.E-3 
BZFRAC(27)=BZFRAC(26)+5.E-3 

DO 71 IZ=28,NZ 

71 BZFRAC(IZ)=BZFRAC(IZ-l)+10.E-3 
DO 72 IZ=1,NZ 

72 BZFRAC(IZ)=BZFRAC(IZ)/GTL 
CALL D0MAIN(1,1,NX,1,NY,1,NZ) 

DO 61 IX=1,NXP1 

DO 62 IY=1,NYP1 
ZL(IX,IY,1)=0.0 
62 ZH(IX,IY,1)=GTL 
DO 63 IZ=1,NZP1 
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YN(IX,IZ.1)=GTH 
63 YS(IX,IZ.1)=0.0 

C YS(IX.13,1) SHOULD COME AFTER 
DO 662 IZ=5,25 
CBL DO 662 IZ=16,25 

662 YS(IX,IZ,1)=(BZFRAC(IZ-1)-BZFRAC(3))xGTABXGTL 
CBL DO 663 IZ=13,15 

CBL GZ12=( BZFRACC IZ-1) -BZFRACC 1 1 ) )XGTL- . 5E-3 
CBL663YS(IX.IZ,1)=SQRTCYS(IX, 16 , 1) XGZ12X2 . -GZ12XX2) 

DO 664 IZ=26.NZ 
664 YS(IX,IZa)=Y$(IX,25,l) 

61 CONTINUE 

STORSACIFIXC LOW) )=. TRUE. 

STORSA(IFIX(HIGH))=.TRUE. 

STORSAC I FIX( SOUTH ))=. TRUE. 

STORWD(IFIX(SOUTH))=.TRUE. 

STORP=.TRUE. 

PRTBFC=.TRUE. 

CDAR DARCY=1.E10 

C 

C— GROUP 7. BLOCKAGE* BLOCK< . F . >, IPLANE, IPWRIT 
C XSET CONSTANT POROSITIES OVER SUB-DOMAINS USING* 

C CALL CONPORCIR, TYPE, VALUE, IXF,IXL,IYF,IYL,IZF,IZL), WHERE* 

C IR=RUN SECTION NUMBER, E.G. 1 FOR RUNl SECTION; 'TYPE'^ EAST, 

C WEST, NORTH, SOUTH, HIGH, LOW & CELL. ' VALUE* =WANTED POROSITY 

C OVER REGION IXF, . . .IZL. 

C ^DIMENSION ARRAYS PECNX, NY, NZ) , PNC NX, NY, NZ) , PHC NX, NY, NZ) , & 

C PC(NX,NY,NZ) ABOVE. 

C XFOR FULLY-BLOCKED CELLS- (IE. *VALUE*= 0.0) USER NEED SET ONLY 
C THE 'CELL* POROSITY (TO ZERO), AS CELL-FACE AREAS ARE THEN 

C AUTOMATICALLY ZEROED. 

C XFOR SATELLITE PRINTOUT OF ALL POROSITIES IN DOMAIN, *IPLANE*= 
C XPLANE YPLANE OR ZPLANE, FOR DESIRED CROSS-SECTION DIRECTION. 
C XFOR EACH 'TYPE* A MAXIMUM OF 10 CALLS TO CONPOR IS ALLOWED, 

C BUT IF REQUIREMENTS EXCEED THIS PROVISION SET BLOCK=.T. & 

C IPWRIT=-1, AND SET POROSITY ARRAYS EXPLICITLY HERE AS WANTED. 
C IN THIS CASE, THE USER MUST SET A L L ELEMENTS OF 

C ARRAYS PE, PN, PH, PC (MANY MAY BE 0.0 OR 1.0). HE MAY USE* 

C CALL CRCPARRAY, VALUE, IXF, IXL , lYF, lYL , IZF, IZL , NX, NY, NZ) 

C ANY NUMBER OF TIMES, TO SET *PARRAY* (= PE, ETC.) TO 
C 'VALUE* OVER RANGE IXF TO IXL, lYF TO lYL, IZF TO IZL. 

C XCONPOR MUST NOT BE USED IN CONJUNCTION WITH EXPLICIT 
C SETTINGS OF THE ARRAYS (INCLUDING SETTINGS VIA CR) . 

C 

C— GROUP 8. DEPENDENT VARIABLES TO BE SOLVED FOR OR STORED * 

C SOLVARC 1-25 )<25X.F.>,ST0VAR( 1-25 )<25X.F.>,C0NC1( 1-4 )<4X.T.> 

C USE FOLLOWING NAMED INTEGERS FOR ARRAY ELEMENTS 1-20: 

C P1,PP,U1,U2,V1,V2,W1,W2,M1,M2,RS,KE,EP,H1,H2,H3,C1,C2,C3,C4. 

S0LVAR(P1)=.TRUE. 

SOLVAR(PP)=.TRUE. 

S0LVAR(V1)=.TRUE. 

S0LVAR(W1)=.TRUE. 

S0LVAR(H1)=.TRUE. 

CT SOLVAR(KE)=.TRUE. 

CT SOLVAR(EP)=.TRUE. 

ST0VAR(V2)=.TRUE. 

ST0VAR(W2)=.TRUE. 

ST0VAR(C1)=.TRUE. 

ST0VAR(C2)=.TRUE. 

ST0VAR(C3)=.TRUE. 

C 

C— GROUP 9. VARIABLE LABELS * 

C TITLE(1-25X2HP1,2HPP,2HU1,2HU2,2HV1,2HV2,2HW1,2HW2,2HR1, 

C 2HR2,2HRS,2HKE,2HEP,2HH1,2HH2,2HH3,2HC1,2HC2, 

C 2HC3,2HC4,2HRX,2HRY,2HRZ, 2x4Hxxxx> 

TITLE(C1)=TITC1 

TITLE(C2)=TITC2 

TITLE(C3)=TITC3 

TITLE(PP)=TITPP 

C 

C— GROUP 10 PROPERTIES: 

C IRH01<1>,IRH02<1>,RH0K1 .0>,RH02<1 .0>, 
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ARH01<1.0>,BRH0K1 .0>,CRH01<1 .0> 
IEMlil<l>» EMUKl . 0>, EMULAM<1 . E-10> 
IHSAT,H1SAT,H2SAT,PSATEX<1 .0> 

SIGMA(1-25X1.0,2.0,1.,1.E10,1.»1.E10,1 

4X1. 0,1. 314, 1.0,1. El 0,1 0X1. 0> 

IRH01=-1 

PT0T=55.E5 

T0T=555.55 

RAIR=287. 

GAMA=1 .35 

CP=RAIR/(1-1/GAMA) 

TW=323. 

HWALL=THXCP 

HT0T=CPXT0T 

RHT0T=PT0T/T0T/RAIR 

L0GIC(87)=.TRUE. 

ARH01=RHT0T/PT0TXX( 1/GAMA) 

BRH01=1 ./GAMA 
TURBULENT OR LAMINAR 
IEMU1=-1 
IEMU1=1 
JEMU1=IEMU1 
EMU1=1 .E-5 
EMULAM=EMU1 
GEMU1=EMU1 
GPR=.7 

SIGMA(24)=GPR 

SIGMA(14)=GPR 



I.EIO, 



— GROUP 11 INTER-PHASE TRANSFER PROCESSES : 

ICFIP,CFIPS, IMDOT,CMDOT,CA1K1 . E6>,CA2K1 . E6> 



— GROUP 12 SPECIAL SOURCES t 

ISPCS0( 1-25) , AGRAVX, AGRAVY, AGRAVZ, ABUOY, HREF 



— GROUP 13 INITIAL FIELDS : 

FI INIT( 1-25X25X1. E-10> 

MACH NO. OF FREE STREAM 
GMACH=3.2 

A=1+(GAMA-1)/2XGMACHXX2 

TE=TOT/A 

RHE = RHTOT/AXX( 1/ ( GAMA-1 ) ) 

PSTAT=PTOT/AXX( GAMA/ (GAMA-1 )) 

RH01=ARH01XPSTATXXBRH01 

SONIC=SQRT(GAMAXRAIRXTE) 

WIN=SONICXGMACH 

RKEIN=0.01XWINXX2 

EPIN=0.16XRKEINXX1 .5/GTH/2. 

FIINIT(W1)=WIN 

FIINIT(P1)=PSTAT 

FIINIT(H1)=HT0T 

FIINIT(KE)=RKEIN 

FIINIT(EP)=EPIN 



— GROUP 14 BOUNDARY/INTERNAL CONDITIONS i 

ILOOP1,ILOOPN,XCYCLE<.F.>,PBAR,REGION(1-10X10X.T.> 

XN.B. ALL 10 REGIONS ARE DEFAULTED .TRUE.. THE USER SHOULD 
SET REGION(I)=. FALSE. FOR UNUSED REGIONS 'I*. 

DO 14 1=1,10 
14 REGION(I)=. FALSE. 



— GROUP 15 TO 24; REGIONS 1 TO 10 
' ONLY THOSE REGIONS ARE ACTIVE WHICH ARE SPECIFIED BY THE 
USER, PREFERABLY BY WAY OFt- 

CALL PLACE(IREGN,TYPE,IXF,IXL,IYF,IYL,IZF,IZL) & 

CALL COVAL(IREGN,VARBLE,COEFF, VALUE) 

CALL PLACE(l,LOW,l,NX,l,NY,l,l) 

C CALL C0VAL(1,M1,FIXFLU,WINXRHE) 

CDAR CALL COVAL ( 1 , Ml , 1 . E-20 , 1 . E+20XWINXRHE) 

GCM=2XGAMA/WI N/ ( GAMA-1 ) 

GVM=PTOTXRHE/RHTOT 
CALL COVAL (1, Ml, GCM,GVM) 
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C 

C 

C 

C 



C 

C 



CT 
CT 
C — 
C — 
C 
C 



CALL COVAL(l,Wl,ONLYMS,WIN) 

CALL C0VAL(1,H1,0NLYMS,HT0T) 

CALL COVAL(l,KE,ONLYMS,RKEIN) 

CALL COVAL(l»EP,ONLYMS,EPIN) 

CALL PLACE(2,HIGH,1,NX,1,NY,NZ,NZ) 

CALL C0VAL(2,M1/FIXVAL,PSTAT*0.) 

CALL COVAL(2,M1,1000*WINkRHE/PSTAT,PSTAT) 
CALL C0VAL(2,H1,0NLYMS,HT0T) 

WALL ALONG THE VANE IZ(11,NZ) 

GCM=EMUl/( . 5XBYFRAC( 1 )XGTH) 
DY1=BYFRAC(1)*GTH 
G0EFF=EMU1/(0.5 xdYI) 
G0EFH=EMU1/(0.5XDY1XSIGMA(2A)) 

CALL PLACE(3,S0UTH,1,NX, 1, 1,A,NZ) 

CALL C0VAL(3,W1,G0EFF,0. ) 

CALL C0VAL(3,H1,G0EFH,HWALL) 

CALL COVAL(3,W1,WALL,0.) 

CALL C0VAL(3,H1,WALL,HWALL) 

CALL C0VAL(3,KE,WALL,0. ) 

CALL COVAL(3,EP,WALL,0.) 



GROUP 25 GROUND STATION i 
GROSTA<.F.>,NAMLST<.F.> 

xnamlst activates namelist in ground. 

GROSTA=.TRUE. 



C 

C GROUP 26 SOLUTION TYPE AND RELATED PARAMETERS 

C WHOL EP< . F . > , SUBPST< . F . > , DONACC< . F . > 

WHOLEP=.TRUE. 



: 



c 

C— GROUP 27 SWEEP AND ITERATION NUMBERS i 
C FSWEEP<1>,LSWEEP<1>,LITHYD<1>AITC<1>,LITKE<1>,LITH<1>. 

C LITER(1-25X9X1,-1,15X1> 

C IVELF<1>,NVEL<1>,IVELL<10000>, 

C IKEF<1>.NKE<1>.IKEL<10000>, 

C IENTF<1>,NENT<1>, IENTL<10000>, 

C ICNCF<1>,NCNC<1>, ICNCL<10000>, 

C IRH01F<1>,NRH01<1>,IRH01L<10000>, 

C IRH02F<1>,NRH02<1>,IRH02L<10000> 

LSWEEP=^00 
GSWP=LSWEEP 
CR FSWEEP=200 

LITER(PP)=20 
LITER(V1)=5 
LITER(W1)=5 
C LITHYD=2 

C 

GROUP 28 TERMINATION CRITERIA : 

C ENDIT ( 1-25 X9X1 .E-10, 0.5,15X1 .E-10> 

ENDIT(l)=l.E-5 

C 

C— GROUP 29 RELAXATION : 

C RLXP<1.>,RLXPXY<1.>,RLXPZ<1.>,RLXRH0<1.>,RLXMDT<1 .>, 

C DTFALSC 3-25X23x1. El 0> 

DTFALS(Wl)=l.E-5 

DTFALS(Vl)=l.E-5 

RLXP=.2 

C 

C-— GROUP 30 LIMITS i 

C VELMAX<1 .E10>,VELMIN<-1 .E10>,RHOMAX<1 .E10>,RHOMIN<1 .E-10>, 
C TKEMAX<1 . E10>, TKEMIN<1 . E-10>, EMUMAX<1 . E10>, EMUMIN<1 . E-10>, 
C EPSMAX<1 . E10>, EPSMIN<1 . E-10>, AMDTMX<1 . E10>, AMDTMN<-1 . E10> 

C— -T 

C— - GROUP 31 SLOWING DEVICES : SL0RH0<1 . >, SL0EMU<1 . > 

SL0RH0=.2 

C 

C— GROUP 32 PRINT-OUT OF VARIABLES : 

C PRINTC1-25X.T. , . F . , 23x . T . >, SUBWGR< . F. > 

PRINTCC1)=.TRUE. 

PRINT(C2)=.TRUE. 

PRINT(C3)=.TRUE. 

PRINT(PP)=.TRUE. 
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C VAN05050 

r GROUP 33 MONITOR PRINT-OUT i VAN05060 

C IXMON<l>,IYMON<l>,IZMON<l>,NPRMON<l>,NPRMNT<l> VAN05070 

NPRM0N=5 VAN05080 

IYM0N=2 VAN05090 

IZM0N=12 VAN05100 

VAN05110 

r— group 34 FIELD PRINT-OUT CONTROL « VAN05120 

C NPRINT<100>,NTPRIN<100>,NXPRIN<1>,NYPRIN<1>,NZPRIN<1>, VAN05130 

C IZPRF<1>,ISTPRF<1>,IZPRL<10000>,ISTPRL<10000> VAN05140 

C NUMCLS<10>,KOUTPT VAN05150 

NPRINT=LSWEEP VAN05160 

C VAN05170 

Q group 35 TABLE CONTROL : VAN05180 

C TABLES<.F.>,NTABLE,NTABVR,LINTAB,NPRTAB,NM0N, VAN05190 

C ITAB(1-8),MTABVR(1-8) VANOSZOO - 

C ^ VAN05210 

C GROUP 36-38 ARE NOT DOCUMENTED IN THE INSTRUCTION VAN05220 

C MANUAL AND ARE INTENDED FOR MAINTENANCE PURPOSES ONLY VAN05230 

C— group 36 DEBUG PRINT-OUT SLAB AND TIME-STEP « VAN05240 

C IZPRK1>, IZPR2<1>,ISTPRK1>,ISTPR2<1> VAN05250 

C VAN05260 

C GROUP 37 DEBUG SWEEP AND SUBROUTINES s VAN05270 

C KEMU,KMAIN,KINDEX,KGEOM,KINPUT,KSODAT,KCOMPF,KSORCE, VAN05280 

C KS0LV1,KS0LV2,KS0LV3,KC0MPP,KADJST,KFLUX,KSHIFT,KDIF, VAN05290 

C KCOMPU,KCOMPV,KCOMPW,KCOMPR,KWALL,KDBRHO<-l>,KDBEXP,KDBMDT VAN05300 

C KDBGEN VAN05310 

C VAN05320 

C— GROUP 38 MONITOR, TEST, AND FLAG : VAN05330 

C M0NITR<.F.>,FLAG<.F.>,TEST<.T.>,KFLAG<1> VAN05340 

C END OF MAINTENANCE-ONLY SECTION VAN05350 

C VAN05360 

C GROUP 39 ERROR AND RESIDUAL PRINT-OUT : VAN05370 

C IERRP<1000>, RESREFCl, 3-24X25X1. >,RESMAP<.F.>, VAN05380 

C RESID(1-25X2X.F. ,23X.T.>,K0UTPT VAN05390 

RESREF(1)=WINXRHE VAN05400 

RESREF(7)=WINXRESREF(1) VAN05410 

RESREF(5)=WINXRESREF(1)X0.1 VAN05420 

RESREF(Hl)=HTOTXRESREF(l) VAN05430 

RESREF(KE)=RKEINXRESREF(1) VAN05440 

RESREF(EP)=EPINXRESREF(1) VAN05450 

IERRP=LSWEEP/10 VAN05460 

KOUTPT=LSWEEP/10 VAN05470 

C VAN05480 

C GROUP 40 SPECIAL DATA : LOGIC( 1 . . 10 ) , INTGR( 1 . . 10 ) , RE( 21 . . 30) , VAN05490 

C NLSP<1>,NISP<1>,NRSP<1>,SPDATA<.F.>,LSPDA(1),ISPDA(1),RSPDA<1) VAN05500 

C USE FIRST 10 ELEMENTS OF ARRAYS LOGIC & INTGR AND 21ST VAN05510 

C TO 30TH OF ARRAY RE FOR TRANSFERRING SPECIAL DATA FROM VAN05520 

C SATELLITE TO GROUND, BUT IF REQUIREMENTS EXCEED THIS VAN05530 

C PROVISION SET SPDATA = .T., AND DIMENSION ARRAYS LSPDA, VAN05540 

C ISPDA, RSPDA ABOVE AND IN GROUND AS NEEDED, AND SET HERE VAN05550 

C VAN05560 

C GROUP 42 RESTARTS AND DUMPS s SAVEM< . F . >, RESTRT< . F . >, KINPUT VAN05570 

SAVEM=.TRUE. VAN05580 

BFPLOT=.TRUE. VAN05590 

C RESTRT=.TRUE. VAN05600 

C •- VAN05610 

C VAN05620 

C GROUP 43 GRAFFIC « VAN05630 

C GRAPHS< . F. >, ORTHOG< . T . >, ANTSYM, NPRT<1>, ITITL<5X4HXXXX> VAN05640 

C— FOR A GRAFFIC RUN, DIMENSION PHIl 8 PHI2 AS FOLLOWS i VAN05650 

C PHIKNXXNYXNZXNM) VAN05660 

C PHI2((NX+2)X(NY+2)X(NZ+2)x<NM+IBLK)) , WHERE VAN05670 

C NM=NO. OF VARIABLES STORED + DENSITY(-IES) VAN05680 

C IBLK=0 IF BL0CK=.FALSE.,=4 IF A 3D RUN, VAN05690 

C =3 IF A 2D.YZ RUN. VAN05700 

C VAN05710 

IF(IRUN.EQ.l) GO TO 900 VAN05720 

900 CONTINUE VAN05730 

C ALL RUNS VAN05740- 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 3 ENDS. VAN05750 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 4 STARTS i VAN05760 
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WRITE GENERAL DATA ON TO THE GUSIEl.DTA TAPE, ETC... 

IF(SPDATA) CALL WRTSPCC LSPDA, NLSP , ISPDA, NISP , RSPDA, NRSP ) 
IF(BLOCK) CALL WRTPORC PE, PN, PH, PC, NX, NY, NZ, IPLANE) 
IF(BFC) CALL WRTBFCC lA , NBFC, XE, XW, YN, YS , ZH, ZL , 

&ND, NX+1 , NY+1 , NZ+1 , NZ, PRTBFC) 

OLD PRACTICES RETAINED FOR REFERENCEi 
IF(SPDATA) CALL SPCDAT(IRUN) 

IF(BLOCK) CALL PORDAT(IRUN) 

IF(GRAPHS) CALL SORT(IRUN) 

IF(RESTRT) GO TO 902 
DO 901 INDVAR=1,25 

IF(IFIX(FIINIT(INDVAR)+0.1).NE. 10101) GO TO 901 
CALL FLDDAT(IRUN) 

GO TO 902 

901 CONTINUE 

902 CALL DATAIO(WRT,10) 

IF(MONITR) CALL DATAIOC WRT , -6 ) 

999 CONTINUE 
STOP 
END 

CXXX IGEN=1 SO BFCXYZ NOT REQUIRED. 

CXXX COMMENT OUT BOTH VERSIONS. 

C 

SUBROUTINE BFCXYZ ( NXPl , NYPl , NZPl ) 

RETURN 

END 



VAN05770 

VAN05780 

VAN05790 

VAN05800 

VAN05810 

VAN05820 

VAN05830 

VAN058<+0 

VAN05850 

VAN05860 

VAN05870 

VAN05880 

VAN05890 

VAN05900 

VAN05910 

VAN059Z0 

VAN05930 

VAN059A0 

VAN05950 

VAN05960 

VAN05970 

VAN05980 

VAN05990 

VAN06000 

VAN06010 

VAN06020 

VAN06030 



56 



Appendix B 
Ground Listing 
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FILEi vantgrd fortran A1 



C$DIRECTIVE^)(MAIN AMI LEITNER 

C LECGRD LAST GEO. NZ=27 NY=18 LAMINAR FLOW 
C XFILE NAME: MODBFCGD.FTN 

C ^aNCLUDE DED SUBROUTINES: THE MODELS OF MAIN, GROUND & STRIDE. 
C ^(DOCUMENTATION: PHOENICS INSTRUCTION MANUAL (SPRING 1983) 

C WITH BODY-FITTED COORDINATES INSTRUCTION SUPPLEMENT 

C (SUMMER 198A). 

C ^SATELLITE FILE NAME: MODSTL.FTN 

C0MM0N/ISHIFT/III(57),NFMAX 

C SET F-ARRAY DIMENSION AS NEEDED, & SET NFMAX ACCORDINGLY. 

C FOR BFC»S ALSO SET Fl-ARRAY DIMENSION AS NEEDED ,AND SET 
C NFIMAX ACCORDINGLY. 

COMMON/FOB/FK 10000) 

COMMON/NFOB/NFIMAX 
COMMON F(25000) 

NFMAX=25000 

NF1MAX=10000 

CALL MAINl 

STOP 

END 

C$DIRECTIVE5()(GROUND 

SUBROUTINE GROUND( IRN, ICHAP , ISTP , ISWP, IZED,INDVAR) 

INCLUDE (CMNGUS) 

INCLUDE (GUSSEQ) 

C INCLUDE NMLIST 

LOGICAL BFC 

EQUIVALENCE ( LOGIC( 20 ) , BFC ) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 1 STARTS: 

C 

C+++++MEANING OF SUBROUTINE ARGUMENTS: 

C IRN=RUN NUMBER; ICHAP=CHAPTER CALLED; ISTP=TIME STEP; 

C ISWP=SOLUTION SWEEP; IZED=Z-SLAB; INDVAR: SEE CHAPTERS BELOW. 

C+++++USER-INTRODUCED VARIABLES & ARRAYS: 

C TO AVOID CONFLICT WITH VARIABLE NAMES USED IN COMMON, ALL 
C VARIABLES INTRODUCED BY THE USER SHOULD HAVE NAMES STARTING 
C WITH »G» IF REAL, IF INTEGER, AND »G» OR IF LOGICAL. 

C THUS GDZ(IZ) MIGHT BE A Z-INTERVAL ARRAY; 

C GW1(IY,IX) A 2-D ARRAY FOR AXIAL VELOCITY; ETC. 

C USER-GENERATED SUBROUTINES SHOULD BE NAMED CORRESPONDINGLY, EG 

C SUBROUTINE GVISC(GTEMP, GCNC, GVSC) , FOR COMPUTING VISCOSITY 
C FROM CONCENTRATION & TEMPERATURE. 

C+++++GROUND-TO-EARTH CONNECTING SUBROUTINES: 

C XUSE GET(NAME,GARRAY,NY,NX) TO PUT VALUES OF VARIABLE NAMED 
C 'NAME* INTO ARRAY 'GARRAY* DIMENSIONED GARRAY( NY, NX) . 

C XUSE SET(NAME,IXF,IXL,IYF,IYL,GARRAY,NY,NX) TO SET VARIABLE 
C 'NAME* TO GARRAY(IY,IX) OVER THE REGION: IXF-IXL 8 lYF-IYL . 

C PRNSLB(NAME) TO PRINT VARIABLE 'NAME* OVER X-Y PLANE. 

C XUSE ADD(NAME,IXF, IXL , lYF, lYL , TYPE, CM, VM, CVAR, WAR, NY, NX) 

C TO ADD SOURCE TO VARIABLE NAMED 'NAME' (SEE CHAPTER 5). 

C ^USE READIZ(IZED) IN CHAPTERS 1, 2, 8, & 9 TO ACCESS PI,.. 

C & VOL,...AHDZ. (SEE FOOTNOTE TO LEGALITY TABLE) 

C 5(USE GET1D(NAME,GARRAY,NDIM) TO PUT VARIABLE NAMED 'NAME* 

C ONE-D ARRAY 'GARRAY' DIMENSIONED NDIM, THUS: 

C CALL GET1D(NAME,GNX,NX) FOR XG,...DXG 8 DIMENSION GNX(NX); 

C CALL GET1D(NAME,GNY,NY) FOR YG,...RV 8 DIMENSION GNY(NY); 

C CALL GET1D(NAME,GNZ,NZ) FOR ZG,...WGRID 8 DIMENSION GNZ(NZ). 

C+++++LEGALITY TABLE FOR USE OF EARTH-CONNECTING SUBROUTINES: 



DM 



IN 



VANOOOlO 

VAN00020 

VAN00030 

VANOOOAO 

VAN00050 

VAN00060 

VAN00070 

VAN00080 

VAN00090 

VANOOlOO 

VANOOllO 

VAN00120 

VAN00130 

VAN00140 

VAN00150 

VAN00160* 

VAN00170 

VAN00180 

VAN00190 

VAN00200 

VAN00210 

VAN00220 

VAN00230 

VAN00240 

VAN00250 

VAN00260 

VAN00270 

VAN00280 

VAN00290 

VAN00300 

VAN00310 

VAN00320 

VAN00330 

VAN00340 

VAN00350 

VAN00360 

VAN00370 

VAN00380 

VAN00390 

VAN00400 

VANOOAIO 

VAN00A20 

VAN00430 

VANOOAAO 

VAN00A50 

VAN00A60 

VAN00A70 

VAN00480 

VAN00490 

VAN00500 

VAN00510 

VAN00520 

VAN00530 

VANOOS^^O 

VAN00550 

VAN00560 

VAN00570 



c 

c 

c 


ENTRIES IN TABLE GIVE 
USED FOR VARIABLES IN 
STRIDE IS REGARDED AS 


CHAPTERS 
LEFT-HAND 
BEING IN 


IN WHICH 
1 COLUMN. 
CHAPTER 


SUBROUTINES 

(SUBROUTINE 

3) 


CAN 1 


BE 


VAN00580 

VAN00590 

VAN00600 

VAN00610 

VAN00620 

VAN00630 

VAN00640 

VAN00650 


U 

c 

c 

p 


: VARIABLE: : GET 8 

: :: PRNSLB 


: SET 

: 


: ADD 
: 


: 


READIZ : 


GETID 




L 

c 


:P1 - RZ : : ALL 


: 6 8 7 


: 5 


. 


1,2, 8, 9: 


NONE 


: 


c 


:P10 - RZH: :3-7, 10-16 


: 3 


: NONE 


: 


NONE : 


NONE 


: 


VAN00660 


c 


:VOL -AHDZ: : ALL 


: 3 


: NONE 


: 


1,2, 8, 9: 


NONE 


: 


VAN00670 


c 


:D1DP :: NONE 


: 10 


: NONE 


: 


NONE : 


NONE 


: 


VAN00680 


c 


|D2DP :: NONE 


: 11 


: NONE 


: 


NONE : 


NONE 


: 


VAN00690 


c 


:MU1,MU1H : : 5,13-16 


: 12 


: NONE 


: 


NONE : 


NONE 


: 


VAN00700 


c 


:EXCO(L,H):: NONE 


: 13 


: NONE 


: 


NONE : 


NONE 


: 


VAN00710 


c 


:CFP :: 5 


: 14 


: NONE 


: 


NONE : 


NONE 


: 


VAN00720 
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FILE: VANTGRD FORTRAN A1 



c iMDT :: 5 : 15 « NONE : NONE : NONE : 

c :HST1,HST2:: 5 & 15 : 16 : NONE : NONE : NONE : 

c :XG -WGRID:: NONE i NONE : NONE t NONE : ALL t 

C 

C NOTES ON ABOVE TABLE: 

C XIN CHAPTERS 1, 2» & 9 VARIABLES PI... DM & GEOMETRY 

C VOL...AHDZ CAN BE ACCESSED BUT ONLY IN CONJUNCTION WITH 
C USE OF READIZ, THUS: 

C DO 1 IZED=1,NZ 

C CALL READIZ(IZED) 

C 1 CALL GET(... AS REQUIRED..) 

C XGEOMETRY ACCESSED BY READIZ IS THAT AT INITIAL TIME. 

C XDIDP & D2DP ONLY ACCESSIBLE IN UNSTEADY FLOWS.' 

C+++++GROUND SERVICE SUBROUTINES: 

C XUSE C0NTUR(NAME,IPLANE,IL0C,NINT,I1,I2, J1,J2,GARRAY,NDIM) FOR 
C LINE-PRINTER PLOTS OF CONTOURS. 'NAME' = U1,...CA; 

C 'IPLANE'= XPLANE, YPLANE, OR ZPLANE; ILOC SETS IX, lY, OR 

C IZ LOCATION OF IPLANE; II, 12, Jl, & J2 SET FIRST & LAST 

C CELLS IN HORIZ. & VERT. ON PLOT; GARRAY IS 1-D WORKING ARRAY 

C OF DIMENSION NXXNY, NXXNZ, OR NYXNZ DICTATED BY IPLANE; & 

C NDIM SETS VALUE OF DIMENSION OF GARRAY. 

C XUSE FLD2DA(TITLE,GARRAY,NY,NX) TO PRINT ANY ARRAY DIMENSIONED 
C GARRAY(NY,NX); SET 'TITLE' TO REQUIRED NAME ( A HOLLERITH 

C CHARACTERS ONLY). 

C XUSE FLD3DA(TITLE, GARRAY, NX, NY, NZ, IPLANE, ILOC) TO PRINT ANY 
C ARRAY DIMENSIONED GARRAYC NX, NY, NZ) IN PLANE SPECIFIED BY 

C 'IPLANE' & 'ILOC AS FOR CONTUR ABOVE; SET 'TITLE' AS FOR 

C FLD2DA. 

C VARIABLE NAMES FOR USE IN GROUND: 

COMMON/TYPE/CELL , EAST, WEST, NORTH, SOUTH, HIGH, LOW, VOLUME, WALL 
C0MM0N/VAR/P1,PP,U1,U2,V1,V2,W1,W2,R1,R2,RS, 
&KE,EP,H1,H2,H3,C1,C2,C3,CA,RX,RY,RZ,S1,S2 
COMMON/VAROLD/PlO,PPO,UlO,U2O,V10,V2O,WlO,W2O,R10,R20,RS0, 
&KEO,EP0,H10,H20,H30,C10,C2O,C30,CAO,RX0,RYO,RZO,S10,S20 
C0MM0N/VARL0W/P1L,PPL,U1L,U2L,V1L,V2L,W1L,W2L,R1L,R2L,RSL, 
&KEL,EPL,H1L,H2L,H3L,C1L,C2L,C3L,CAL,RXL,RYL,RZL,S1L,S2L 
COMMON/VARHI/P1H,PPH,U1H,U2H,V1H,V2H,W1H,W2H,R1H,R2H,RSH, 
&KEH,EPH,H1H,H2H,H3H,C1H,C2H,C3H,CAH,RXH,RYH,RZH,S1H,S2H 
COMMON/GMTRY/VOL,VOLO,AEAST,ANORTH,AHIGH,AEDX,ANDY, AHDZ 
COMMON/PROP/Dl, D2,D1DP, D2DP,MU1 , MUILAM, EXCO,CFP, MDT, HSTl, HST2 
COMMON/PRPOLD/DIO, D20 
COMMON/PRPLOW/DIL , D2L , EXCOL 
COMMON/PRPHI/D1H,D2H,MU1H, EXCOH 
COMMON/VARNX/XG,XU, DXU, DXG 
COMMON/VARNY/YG,YV,DYV,DYG,R,RV 
COMMON/ VARNZ/ZG, ZWl , DZW, DZG, WGRID 
COMMON/GDMSCI/XPLANE, YPLANE, ZPLANE, ITNO 
COMMON/GDMSCL/LSLAB,MSLAB,HSLAB,LAMMU 
REAL NORTH, LOW 

INTEGER P1,PP,U1,U2,V1,V2,W1,W2,R1.R2,RS, 

&EP,H1,H2,H3,C1,C2,C3,CA,RX,RY,RZ,S1,S2 
INTEGER P10,PPO,U10,U20,V10,V20,W10,W20,R10,R20,RSO, 

&EP0,H10,H2O,H30,C10,C20,C30,CA0,RX0,RY0,RZO,S10,S20 
INTEGER P1L,PPL,U1L,U2L,V1L,V2L,W1L,W2L,R1L,R2L,RSL, 

&EPL,H1L,H2L,H3L,C1L,C2L,C3L,CAL,RXL,RYL,RZL,S1L,S2L 
INTEGER P1H,PPH,U1H,U2H,V1H,V2H,W1H,W2H,R1H,R2H,RSH, 

&EPH,H1H,H2H,H3H,C1H,C2H,C3H,CAH,RXH,RYH,RZH,S1H,S2H 
INTEGER VOL, VOLO,AEAST,ANORTH,AHIGH,AEDX, ANDY, AHDZ 
INTEGER D1 , DIDP, D2, D2DP, EXCO, CFP, HSTl , HST2 
INTEGER DIO, D20, DIL , D2L, EXCOL , DIH, D2H, EXCOH 
I NT EG ER XG , XU , DXU , DXG , YG , YV , DYV , DYG , R , RV , ZG, ZWl , DZW , 

&DZG, WGRID 

INTEGER XPLANE, YPLANE, ZPLANE 
LOGICAL LSLAB,MSLAB,HSLAB,LAMMU,LSPDA 
EQUIVALENCE ( Ml , R1 ) , (M2, R2) 

C SATLIT-EQUIVALENT IRUN: 

EQUIVALENCE ( IRUN, INTGR( 11 ) ) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 1 ENDS. 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 1 STARTS: 

C ARRAYS ( DIMENSIONED NY, NX ) FOR USE WITH 'ADD': 

DIMENSION CVAR(l,l),VVAR(l,l),CM(l,l),VM(l,l),ZERO(l,l) 
DIMENSION GP(30,1),GH(30,1),GD(30,1),GV(30,1),GW(30,1) 



VAN00730 

VAN007A0 

VAN00750 

VAN00760 

VAN00770 

VAN00780 

VAN00790 

VAN00800 

VAN00810 

VAN00820 

VAN00830 

VAN008A0 

VAN00850 

VAN00860 

VAN00870 

VAN00880 

VAN00890 

VAN00900 

VAN00910 

VAN00920 

VAN00930 

VAN009A0 

VAN00950 

VAN00960 

VAN00970 

VAN00980 

VAN00990 

VANOIOOO 

VANOlOlO 

VAN01020 

VAN01030 

VANOIOAO 

VAN01050 

VAN01060 

VAN01070 

VAN01080 

VAN01090 

VANOllOO 

VANOlllO 

VAN01120 

VAN01130 

VANOllAO 

VAN01150 

VAN01160 

VAN01170 

VAN01180 

VAN01190 

VAN01200 

VAN01210 

VAN01220 

VAN01230 

VAN012A0 

VAN01250 

VAN01260 

VAN01270 

VAN01280 

VAN01290 

VAN01300 

VAN01310 

VAN01320 

VAN01330 

VAN013A0 

VAN01350 

VAN01360 

VAN01370 

VAN01380 

VAN01390 

VANOIAOO 

VANOIAIO 

VAN01A20 

VAN01A30 

VANOIAAO 



59 



FILEi VANTGRD FORTRAN A1 



1 ,GMACH(30,1)»GTEMP(30,1),GVISC(30,1),GWH(30,1),GWM(30,1) VAN01A50 

2 ,GKE(30,1)»GC3(30,1),GYG(30,1),GXX(30,1),GYY(30, 1)»GZZ(30,1) VAN01A60 

C SPECIAL-DATA ARRAYS DIMENSIONED & DIMENSION VALUES SET HERE* VAN01A70 

DIMENSION LSPDA(1),ISPDA(1),RSPDA(1) VAN01A80 

C USER PLACES HIS VARIABLES, ARRAYS, EQUIVALENCES ETC. HERE. VAN01490 

EQUIVALENCE ( RAIR, RE( 21) ) , (GAMA, RE( 22 ) ) , (GSWP, RE( 23) ) , VAN01500 

1(GPR,RE(24)),(GTW,RE(25)),(GEMU1,RE(26)),(JEMU1,INTGR(D) VAN01510 

DATA NLSP,NISP,NRSP/1,1,1/ VAN01520 

DATA CVAR,VVAR,CM,VM,ZERO/5XO.O/ VAN01530 

C USER PLACES HIS DATA STATEMENTS HERE. VAN01540 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 1 ENDS. VAN01550 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 2 STARTS: VAN01560 

C PLEASE DO NOT ALTER, OR RE-SET, ANY OF THE REMAINING VAN01570 

C STATEMENTS OF THIS SECTION. VAN01580 

DATA NUMCH4 / 0 / VAN01590 

IF(SPDATA) VANO16a0' 

aCALL RDSPC(IRN,INTGR(12),LSPDA,NLSP,ISPDA,NISP,RSPDA,NRSP) VAN01610 

CALL GRDUTY(IRN,ICHAP,IZED,INDVAR) VAN01620 

IF(BFC) CALL BFCGRDC IRN, ICHAP, ISWP, IZED, INDVAR) VAN01630 

IF(ICHAP.EQ.-5) GO TO 10 VAN01640 

IF(ICHAP.LE.0.0R.ICHAP.GT.16) RETURN VAN01650 

GO TO (100,200,300,4999,500,600,700,800,900,1000,1100,1200, VAN01660 

81300,1400,1500, 1600), ICHAP VAN01670 

RETURN VAN01680 

4999 NUMCH4= NUMCH4 + 1 VAN01690 

IF (M0D(NUMCH4,2) .EQ.l) GO TO 400 VAN01700 

RETURN VAN01710 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX STANDARD SECTION 2 ENDS. VAN01720 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION 2 STARTS: VAN01730 

C VAN01740 

C CHAPTER 0: MODIFY SATLIT DATA, AT START OF EACH IRN. VAN01750 

C VAN01760 

10 CONTINUE VAN01770 

C IF( .NOT.NAMLST) RETURN VAN01780 

C IF(IRN.EQ.NRUN) DATFIL= . FALSE. VAN01790 

C— READ SATLIT DATA NAMELIST HERE VAN01800 

C CALL WRIT40(40HENTER NAMELIST DATA FOR GROUPS 1 TO 24 ) VAN01810 

C READ(20,G1G24) VAN01820 

C CALL WRIT40(40HENTER NAMELIST DATA FOR GROUPS 25 TO 42 ) VAN01830 

C READ(20,G25G42) VAN01840 

RETURN VAN01850 

C VAN01860 

C CHAPTER 1: CALLED AT THE START OF EACH TIME STEP. VAN01870 

C SET 'DT* HERE WHEN TLAST SET NEGATIVE IN BLOCK DATA. VAN01880 

C 'ATIME + DT* GIVES THE END TIME OF THE CURRENT TIME STEP. VAN01890 

C NOT ACCESSED IF STEADY, OR PARABOLIC. VAN01900 

0 VAN01910 

100 CONTINUE VAN01920 

RETURN VAN01930 

0 VAN01940 

C CHAPTER 2: CALLED AT THE START OF EACH SWEEP. VAN01950 

0 VAN01960 

200 CONTINUE VAN01970 

RETURN VAN01980 

0 VAN01990 

C CHAPTER 3: CALLED AT THE START OF EACH SLAB; VAN02000 

C NOT ACCESSED IF PARABOLIC,' BUT 'STRIDE* IS. VAN02010 

0 VAN02020 

300 CONTINUE VAN02030 

RETURN VAN02040 

0 VAN02050 

C CHAPTER 4: CALLED AT THE START OF EACH RE-CALCULATION OF VAN02060 

C VARIABLES P1,...C4 AT CURRENT SLAB. ITNO= ITERATION NUMBER. VAN02070 

0__L VAN02080 

400 CONTINUE VAN02090 

RETURN VAN02100 

0 VAN02110 

C CHAPTER 5: GROUND CALLED WHEN SOURCE TERM IS COMPUTED. VAN02120 

C INDVAR GIVES DEPENDENT VARIABLE IN QUESTION IE. U1,...C4. VAN02130 

C TO ADD SOURCE TO DEPENDENT VARIABLE CKSAY) FOR IX=IXF,IXL VAN02140 

C AND IY=IYF,IYL INSERT STATEMENT: VAN02150 

C IF(INDVAR.EQ.Cl) VAN02160 
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ooo oooo oooooooooooooooooo 



FILEi vantgrd fortran ai 



&CALL ADD(INDVAR,IXF,IXL,IYF, IYL,TYPE,CM,VM,CVAR,VVAR,NY,NX) 
NOTES ON *ADD'« 

XS0URCE= (CVAR(IY,IX)+AMAX1(0.0,MASFL0))X(VVAR(IY,IX)-PHI), 

WHERE 'PHI' IS IN-CELL VALUE OF VARIABLE IN QUESTION. 

X'MASFL0'= CM(IY,IX)X(VM(IY,IX)-P), 

WHERE 'P' IS THE IN-CELL PRESSURE. 

XFOR INDVAR= Ml, OR =M2, SOURCE ADDED IS 'MASFLO' ONLY, 

EXCEPT FOR ONEPHS=.F. & MASFLO < 0.0 (IE. OUTFLOW) WHEN 
CM(IY,IX) IS MULTIPLIED BY RlxDl (FOR Ml) S, R2XD2 (FOR M2). 

XBOTH 'CVAR' & 'CM' ARE MUTLIPLIED BY CELL-GEOMETRY QUANTITY 
DICTATED BY SETTING OF 'TYPE' (=CELL, EAST AREA, .. VOLUME) . 

XTYPE-SPECIFIED AREAS ARE CALCULATED AS IF BLOCKAGE ABSENT, 

BUT 'VOLUME' WITH ACCOUNT FOR ITS PRESENCE. 

XFOR ALL SOLVED VARIABLES, INCLUDE DING Ml ( & M2 WHEN 0NEPHS=F), 
IF 'CM'> 0.0 CALL 'ADD'; FOR Ml & M2 ALTHOUGH 'CVAR' & 'WAR' 
HAVE NO SIGNIFICANCE THEY MUST BE ENTERED AS ARGUMENTS. 

X'CVAR', 'WAR', 'CM' & 'VM' MUST BE DIMENSIONED NY, NX. 



500 CONTINUE 
RETURN 



CHAPTER 6s CALLED AT THE END OF EACH VARIABLE-RECALCULATION 
CYCLE COMMENCED AT CHAPTER 4. ITNO = ITERATION NUMBER. 



600 CONTINUE 
RETURN 



CHAPTER 7: CALLED AT END OF EACH SLAB-WISE CALCULATION. 



700 CONTINUE 

IF(FLOATdSWP) .LT.GSWP) RETURN 
CALL GET(P1,GP,NY,NX) 

CALL GET(H1,GH,NY,NX) 

CALL GET(D1,GD,NY,NX) 

CALL GET(V1,GV,NY,NX) 

CALL GET(W1,GW,NY,NX) 

CALL GET(KE,GKE,NY,NX) 

C CALL GET1D(YG,GYG,NY) 

CALL GRED1(39,IZED,GYG,NY,NX) 

CALL GRED3(57,IZED,GXX,GYY,G2Z,NY,NX) 

GCP=RAIR/(1. -1/GAMA) 

DO 701 1=1, NY 

GSON=SQRT( GAMAXGP( I , 1 )/GD( 1,1)) 
GAV=SQRT(GV(I,1)XX2+GW(I,1)XX2) 

GMACH(I,l)=GAV/GSON 

701 GTEMP(I,1)=GP(I,1)/GD(I,1)/RAIR 

C 701 GTEMP(I>l)=(GH(I,l)-GW(I,l)xx2/2.-GV(I,l)XX2/2.)/GCP 



CALL SET(C1,1,NX,1,NY,GMACH,NY,NX) 

CALL SET(C2,1,NX,1,NY,GTEMP,NY,NX) 

C CALCULATE DYl CF ST H(CONVECTIVE COEF.) Q TAU TR 

IF(JEMU1 .NE.2) GOTO 702 

C TURBULENT VALUES 

GCF=2./GW(NY,1)xx2XGKE(1,1)/3.33XGD(1,1)/GD(NY,1) 

C7 GCF=GCFXGD( NY, 1 )/GD( 1 , 1 )XGTEMP( NY, 1 )/GTEMP( 1 , 1 )XGP( 1 , 1 )/GP( NY, 1 ) 

GST=GCF/2./GPRXX.666 
GHH=GD(NY,1)XGCPXGW(NY,1)XGST 
GR=GPRXX.333 

GTR=GTEMP( NY, 1 )X( 1 , +GRX( GAMA-1 . )/2 . XGMACH( NY, 1 )XX2) 

C !(!.+( GAMA-1 . )/2.xGMACH(NY,1)xx2) 

GQ=GHHX(GTR-GTW) 

GOTO 703 

C LAMINAR VALUES 



702 CONTINUE 

IF(JEMUl.EQ.-l) GEMUl=GVISC(l,l) 
GQ=GEMU1/GPRX(GH(1,1)-GTWXGCP)/GYG(1,1) 

GR=GPRXX . 5 

GTR=GTEMP(NY, 1 )X(1 . +GRX(GAMA-1 . )/2 . XGMACH( NY, 1 )XX2) 
C 1(1 .+(GAMA-1. )/2.XGMACH(NY,l)XX2) 

GHH=GQ/(GTR-GTW) 

GST=GHH/(GD(NY,1)XGW(NY,1)XGCP) 

GTAU=GEMU1XGW( 1 , 1 )/GYG( 1 , 1 ) 
GCF=GTAUX2./(GD(NY,1)XGW(NY,1)XX2) 



VAN02170 

VAN02130 

VAN02190 

VAN02200 

VAN02210 

VAN02220 

VAN02230 

VAN02240 

VAN02250 

VAN02260 

VAN02270 

VAN02280 

VAN02290 

VAN02300 

VAN02310 

VAN02320 

VAN02330 

VAN02340 

VAN02350 

VAN02360 

VAN02370 

VAN02380 

VAN02390 

VAN02400 

VAN02410 

VAN02420 

VAN02430 

VAN02440 

VAN02450 

VAN02460 

VAN02470 

VAN02480 

VAN02490 

VAN02500 

VAN02510 

VAN02520 

VAN02530 

VAN02540 

VAN02550 

VAN02560 

VAN02570 

VAN02580 

VAN02590 

VAN02600 

VAN02610 

VAN02620 

VAN02630 

VAN02640 

VAN026-50 

VAN02660 

VAN02670 

VAN02680 

VAN02690 

VAN02700 

VAN02710 

VAN02720 

VAN02730 

VAN02740 

VAN02750 

VAN02760 

VAN02770 

VAN02780 

VAN02790 

VAN02800 

VAN02810 

VAN02820 

VAN02830 

VAN02840 

VANO2850 

VAN02860 

VAN02870 

VAN02880 
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FILEi VANTGRD FORTRAN A1 



703 GC3(1,1)=GYG(1,1) 

GC3(2,1)=GCF 

GC3C3, 1)=GST 

GC3(<^,1)=GCF/2./GST 

GC3(5,1)=GHH 

GC3(6,1)=GQ 

GC3(7,1)=GTAU 

GC3(8,1)=GTR 

GC3(9,1)=GTR"GTW 

GC 3 ( 1 0 , 1) = G D ( N Y , 1) XGW ( NY , 1) GZZ ( 1 , 1 ) / G EMU 1 
GC3(lia)=GZZ(l,l) 

GC3(12,1)=GEMU1 

GC3(13,l)=GD(NYa)XGW(NY,l)XGYG(l,l)/GEMUlXSQRT(ABS(GCF/2,)) 
CALL SET(C3,l.NXa,NY,GC3,NY,NX) 

RETURN 



CHAPTER 8: CALLED AT THE END OF EACH SWEEP; 
NOT ACCESSED IF PARABOLIC. 



800 CONTINUE 
RETURN 



CHAPTER 9* CALLED AT THE END OF EACH TIME STEP; 
NOT ACCESSED IF PARABOLIC, 



900 CONTINUE 
RETURN 



CH/*^TER lOi SET PHASE 1- DENSITY HERE WHEN IRH01 = -1 IN DATA. 
SET CURRENT-Z 'SLAB' DENSITY, D1 , IF MSLAB=.T., 

EG. IF(MSLAB) CALL SET( D1 , 1 , NX, 1 , NY, GDI , NY, NX) . 

SET NEXT LARGER-Z 'SLAB' DENSITY, DIH, IF HSLAB=,T. 8 PARAB=F 
EG. IFCHSLAB) CALL SET( DIH, 1 , NX, 1 , NY,GD1H, NY, NX) . 

SET D(LN(D1))/DP (IE. DIDP) FOR UNSTEADY FLOW, 

EG. IF(MSLAB) CALL SET( D1 DP , 1 , NX, 1 , NY, GDIDP, NY, NX) . 



1000 CONTINUE 

IF (MSLAB) GO TO 101 

JP1=P1H 

JH1=H1H 

JD1=D1H 

JW1=W1H 

JV1=V1H 

GO TO 102 

101 JP1=P1 
JH1=H1 
JD1=D1 • 

JW1=W1 

JV1=V1 

102 CALL GET(JP1,GP,NY,NX) 

CALL GET(JH1,GH,NY,NX) 

CALL GET(JW1,GW,NY,NX) 

CALL GET( JV1,GV,NY,NX) 

IF(IZED.EQ.l) GOTO 105 
IF(IZED.EQ.NZ) GOTO 109 

C IZED = 2,NZ"1 

DO 103 IX=1,NX 

DO 103 IY=1,NY 

IFCHSLAB) GOTO 10<+ 

GWA=(GW(IY,IX)+GWM(IY,IX))/2. 

GWM(IY,IX)=GW(IY,IX) 

GOTO 115 

104 GWA=(GW(IY,IX)+GWH(IY,IX))/2. 
GWH(IY,IX)=GW(IY,IX) 

115 GHS=GH(IY,IX)-(GWAXX2+GV(IY,IX)XX2)/2. 
IFCGHS.LE.l .E5) GHS=1.E5 

103 GD(IY,IX)= GP(IY,IX)/(1-1/GAMA)/GHS 
GOTO 113 

C IZED = 1 

105 DO 106 IX=1,NX 
DO 106 IY=1,NY 

GHS=GH(IY,IX)-(GW(IY,IX)XX2+GV(IY,IX)XX2)/2. 



VAN02890 

VAN02900 

VAN02910 

VAN02920 

VAN02930 

VAN02940 

VAN02950 

VAN02960 

VAN02970 

VAN02980 

VAN02990 

VAN03000 

VAN03010 

VAN03020 

VAN03030 

VAN0304.0* 

VAN03050 

VAN03060 

VAN03070 

VAN03080 

VAN03090 

VAN03100 

VAN03110 

VAN03120 

VAN03130 

VAN03140 

VAN03150 

VAN03160 

VAN03170 

VAN03180 

VAN03190 

VAN03200 

VAN03210 

VAN03220 

VAN03230 

VAN03240 

VAN03250 

VAN03260 

VAN03270 

VAN03280 

VAN03290 

VAN03300 

VAN03310 

VAN03320 

VAN03330 

VAN03340 

VAN03350 

VAN03360 

VAN03370 

VAN0338'0 

VAN03390 

VAN03400 

VAN03410 

VAN03420 

VAN03430 

VAN03440 

VAN03450 

VAN03460 

VAN03470 

VAN03480 

VAN03490 

VAN03500 

VAN03510 

VAN03520 

VAN03530 

VAN03540 

VAN03550 

VAN03560 

VAN03570 

VAN0358Cr 

VAN03590 

VAN03600 
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FILEi vantgrd fortran A1 



IF(GHS.LE.1.E5) GHS=1.E5 VAN03610 

GD(IY,IX)= GP(IY,IX)/(1-1/GAMA)/GHS VAN03620 

IF(HSLAB) GOTO 107 VAN03630 

GWM(IY,IX)=GW(IY,IX) VAN03640 

GOTO 106 VAN03650 

107 GWH(IY,IX)=GW(IY,IX) VAN03660 

106 CONTINUE VAN03670 

goto 113 VAN03680 

Q IZED=NZ VAN03690 

109 DO 110 IX=1,NX VAN03700 

DO 110 IY=1,NY VAN03710 

IF(HSLAB) GOTO 111 VAN03720 

GHS=GH(IY,IX)-(GWM(IY,IX)XX2+GV(IY,IX)XX2)/2. VAN03730 

IF(GHS.LE.1.E5) GHS=1.E5 VAN037A0 

GWM(IY,IX)=GW(IY,IX) VAN03750 

GOTO 112 VAN0376l0 

111 GHS=GH(IY,IX)-(GWH(IY,IX)XX2+GV(IY,IX)x*2)/2. VAN03770 

C IF(GHS.LE.l .E5) GHS=1.E5 VAN03780 

GWH(IY,IX)=GW(IY,IX) VAN03790 

112 GD(IY,IX)= GP(IY,IX)/(1-1/GAMA)/GHS VAN03800 

110 CONTINUE VAN03810 

C VAN03820 

113 CONTINUE VAN03830 

CALL SET(JD1,1,NX,1,NY,GD,NY,NX) VAN03890 

RETURN VAN03850 

C VAN03860 

C CHAPTER 11: SET PHASE 2 DENSITY HERE WHEN IRH02=-1 IN DATA. VAN03870 

C SET CURRENT-Z 'SLAB' DENSITY, D2, IF MSLAB=.T., VAN03880 

C EG. IF(MSLAB) CALL SET( D2, 1 , NX, 1 , NY, GD2, NY, NX) . VAN03890 

C SET NEXT LARGER-Z 'SLAB* DENSITY, D2H, IF HSLAB=.T. & PARAB=F VAN03900 

C EG. IF(HSLAB) CALL SET( D2H, 1 , NX, 1 , NY,GD2H, NY, NX) . VAN03910 

C SET D(LN(D2))/DP FOR UNSTEADY FLOW, VAN03920 

C EG. IF(MSLAB) CALL SET( D2DP, 1 , NX, 1 , NY, GD2DP, NY, NX) . VAN03930 

C VAN039A0 

1100 CONTINUE VAN03950 

RETURN VAN03960 

C VAN03970 

C CHAPTER 12: SET PHASE 1 VISCOSITY HERE WHEN IEMU1=-1 IN DATA. VAN03980 

C SET CURRENT-Z 'SLAB* VISCOSITY (MUD, IF MSLAB=.T., VAN03990 

C EG. IF(MSLAB) CALL SET(MU1 , 1 , NX, 1 , NY, GVISC, NY, NX) . VANOAOOO 

C SET NEXT LARGER-Z 'SLAB' VISC. (MUIH), IF HSLAB=.T. & PARAB=F VANOAOlO 

C EG. IF(HSLAB) CALL SET(MU1H, 1 , NX, 1 , NY, GVSCH, NY, NX) . VAN09020 

C VAN0A030 

C CHAPTER ALSO ACCESSED WHEN EMULAM=-1.0 IN DATA, SO THAT THE VANOAOAO 

C LAMINAR VISCOSITY WHICH APPEARS IN WALL FUNCTIONS & IN THE VAN0A050 

C KE-EP TURBULENCE MODEL (IEMU1=2) MAY BE SET NON-CONSTANT. VAN0A060 

C SET CURRENT-Z 'SLAB' VALUE (MUILAM) WHEN LAMMU=.T., VAN09070 

C EG. IF(LAMMU) CALL SETCMUILAM, 1 , NX, 1 , NY,GVSCL , NY, NX) . VAN0A080 

C VAN09090 

1200 CONTINUE VAN09100 

GCP=RAIR/(1. -1/GAMA) VANOAllO 

IF (HSLAB) GOTO 122 VAN09120 

CALL GET(H1,GH,NY,NX) VAN09130 

CALL GET(W1,GW,NY,NX) VAN0<4l90 

CALL GET(V1,GV,NY,NX) VAN09150 

GOTO 123 VAN0<il60 

122 CALL GET(H1H,GH,NY,NX) VAN09170 

CALL GET(W1H,GW,NY,NX) VAN09180 

CALL GET(V1H,GV,NY,NX) VAN0<4l90 

123 CONTINUE VAN09200 

DO 121 IX=1,NX VAN0^210 

DO 121 IY=1,NY VAN04220 

GTMP = (GH(IY,IX)-GW(IY,IX)XX2/2.-GV(IY,IX)?<X2/2.)/GCP VAN09230 

IFCGTMP.lt. 150.) GTMP=150. VAN0A2<i0 

121 GVISC(IY,IX)=1.716E-05X(GTMP/273.)XX0.666 VAN04250 

C121 IF(GVISC(IY,IX) .LE. .8E-5) GVISC( lY, IX)= . 8E-5 VAN09260 

IF (MSLAB) CALL SET(MU1 , 1 , NX, 1 , NY, GVISC, NY, NX) VAN04270 

IF (HSLAB) CALL SETCMUIH, 1 , NX, 1 , NY, GVISC, NY, NX) VAN04280 

IF (LAMMU) CALL SETCMUILAM, 1 , NX, 1 , NY, GVISC, NY, NX) VAN09290 

RETURN VAN0A300 

C VAN0A310 

C CHAPTER 13: SET EXCHANGE COEFFICIENT (E.C.) FOR VARIABLE VAN0A320 
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FILEi VANTGRD FORTRAN A1 



INDVAR WHEN SIGMAC INDVAR ) . 0 IN DATA. 

SET CURRENT-Z 'SLAB* E.C. (EXCO) IF MSLAB=.T., 

EG. IF(MSLAB) CALL SET( EXCO ,1 , NX, 1 , NY, GEXCO , NY, NX) . 

SET NEXT SMALLER-Z »SLAB' E.C. (EXCOL) IF LSLAB=.T., 

EG. IF(LSLAB) CALL SET( EXCOL , 1 , NX, 1 , NY, GEXCOL , NY, NX) . 

SET NEXT LARGER-Z 'SLAB' E.C. (EXCOH) IF HSLAB=.T., 

EG. IF(HSLAB) CALL SET( EXCOH, 1 , NX, 1 , NY, GEXCOH, NY, NX) . 

NOTE: FOR MSLAB, INDVAR=U1 , . . CA ; FOR LSLAB, INDVAR=U1L, . .CAL 
& FOR HSLAB, INDVAR=U1H, . . CAH . IF PARAB=.T. SET MSLAB ONLY. 



1300 CONTINUE 
RETURN 



CHAPTER lA; SET INTER-PHASE FRICTION COEFFICIENT (CFP) HERE 
WHEN ICFIP = -1 IN DATA; ITS UNITS = FORCE / (CELL X RELATIVE 
SPEED OF PHASES). 



lAOO CONTINUE 
RETURN 



CHAPTER 15: SET INTER-PHASE MASS-TRANSFER RATE PER CELL (MDT) 
HERE WHEN IMDOT = -1 IN DATA. 



1500 CONTINUE 
RETURN 



CHAPTER 16: SET HERE PHASE 1 a 2 SATURATION ENTHALPIES 
( HSTl a HST2) WHEN IHSAT = -1 IN DATA. 



1600 CONTINUE 
RETURN 
END 



VAN0A330 

VAN0A3A0 

VAN0A350 

VAN0A360 

VAN0A370 

VAN0A380 

VAN0A390 

VANOAAOO 

VANOAAIO 

VAN0AA20 

VAN0AA30 

VANOAAAO 

VAN0AA50 

VAN0AA60 

VAN0AA70 

VAN0AA80 

VAN0AA90 

VAN0A500 

VAN0A510 

VAN0A520 

VAN0A530 

VAN0A5A0 

VAN0A550 

VAN0A560 

VAN0A570 

VAN0A580 

VAN0A590 

VAN0A600 

VAN0A610 

VAN0A620 

VAN0A630 

VAN0A6A0 
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